Есть ли способ уменьшить сумму 100M float элементов массива в CUDA?

Question

Я новичок в CUDA.Поэтому, пожалуйста, имейте в виду вопросы с тривиальными решениями, если таковые имеются.

Я пытаюсь найти сумму 100M плавающих элементов массива.Из следующего кода видно, что я использовал сокращенное ядро ​​и тягу. Я полагаю, ядро ​​хранит сумму в g_odata[0].Поскольку все элементы одинаковы в g_idata, результат должен быть n*g_idata[1].Но вы могли ясно видеть, что результаты неверны для них обоих.

1. Что я делаю не так?Как я могу достичь своей цели?
2. Каждое ядро ​​сокращения, которое я нашел, предназначено для целочисленного типа данных.Например, настоятельно рекомендуется Оптимизация параллельного сокращения в CUDA. .Есть ли какая-то конкретная причина для этого?

Вот мой код:

    #include <iostream>
    #include <math.h>
    #include <stdlib.h>
    #include <iomanip>
    #include <thrust/reduce.h>
    #include <thrust/execution_policy.h>


    using namespace std;


    __global__ void reduce(float *g_idata, float *g_odata) {

    __shared__ float sdata[256];


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

    sdata[threadIdx.x] = g_idata[i];

    __syncthreads();

    for (int s=1; s < blockDim.x; s *=2)
    {
        int index = 2 * s * threadIdx.x;;

        if (index < blockDim.x)
        {
            sdata[index] += sdata[index + s];
        }
        __syncthreads();
    }


    if (threadIdx.x == 0)
        atomicAdd(g_odata,sdata[0]);
    }




    int main(void){

    unsigned int n=pow(10,8);
    float *g_idata, *g_odata;

    cudaMallocManaged(&g_idata, n*sizeof(float));
    cudaMallocManaged(&g_odata, n*sizeof(float));

    int blockSize = 32;
    int numBlocks = (n + blockSize - 1) / blockSize;

    for(int i=0;i<n;i++){g_idata[i]=6.1;g_odata[i]=0;}


    reduce<<<numBlocks, blockSize>>>(g_idata, g_odata);
    cudaDeviceSynchronize();


    cout << g_odata[0] << "\t" << (float)n*g_idata[1] << "\t"<< (float)n*g_idata[1]-g_odata[0]<<endl;

    g_odata[0]=thrust::reduce(thrust::device, g_idata, g_idata+n);

    cout << g_odata[0] << "\t" << (float)n*g_idata[1] << "\t"<< (float)n*g_idata[1]-g_odata[0]<<endl;



    cudaFree(g_idata);
    cudaFree(g_odata);

    }

Результат:

6.0129e+08  6.1e+08 8.7097e+06
6.09986e+08 6.1e+08 13824

Я использую CUDA 10. nvcc --version:

nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2018 NVIDIA Corporation
Built on Sat_Aug_25_21:08:01_CDT_2018
Cuda compilation tools, release 10.0, V10.0.130

Подробная информация о моем GPU DeviceQuery:

./deviceQuery Starting...

 CUDA Device Query (Runtime API) version (CUDART static linking)

Detected 1 CUDA Capable device(s)

Device 0: "GeForce GTX 750"
  CUDA Driver Version / Runtime Version          10.0 / 10.0
  CUDA Capability Major/Minor version number:    5.0
  Total amount of global memory:                 1999 MBytes (2096168960 bytes)
  ( 4) Multiprocessors, (128) CUDA Cores/MP:     512 CUDA Cores
  GPU Max Clock rate:                            1110 MHz (1.11 GHz)
  Memory Clock rate:                             2505 Mhz
  Memory Bus Width:                              128-bit
  L2 Cache Size:                                 2097152 bytes
  Maximum Texture Dimension Size (x,y,z)         1D=(65536), 2D=(65536, 65536), 3D=(4096, 4096, 4096)
  Maximum Layered 1D Texture Size, (num) layers  1D=(16384), 2048 layers
  Maximum Layered 2D Texture Size, (num) layers  2D=(16384, 16384), 2048 layers
  Total amount of constant memory:               65536 bytes
  Total amount of shared memory per block:       49152 bytes
  Total number of registers available per block: 65536
  Warp size:                                     32
  Maximum number of threads per multiprocessor:  2048
  Maximum number of threads per block:           1024
  Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
  Max dimension size of a grid size    (x,y,z): (2147483647, 65535, 65535)
  Maximum memory pitch:                          2147483647 bytes
  Texture alignment:                             512 bytes
  Concurrent copy and kernel execution:          Yes with 1 copy engine(s)
  Run time limit on kernels:                     Yes
  Integrated GPU sharing Host Memory:            No
  Support host page-locked memory mapping:       Yes
  Alignment requirement for Surfaces:            Yes
  Device has ECC support:                        Disabled
  Device supports Unified Addressing (UVA):      Yes
  Device supports Compute Preemption:            No
  Supports Cooperative Kernel Launch:            No
  Supports MultiDevice Co-op Kernel Launch:      No
  Device PCI Domain ID / Bus ID / location ID:   0 / 1 / 0
  Compute Mode:
     < Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >

deviceQuery, CUDA Driver = CUDART, CUDA Driver Version = 10.0, CUDA Runtime Version = 10.0, NumDevs = 1
Result = PASS

Заранее спасибо.

Answer 1

Я думаю, что причина, по которой вы запутались в результатах, заключается в недостаточном понимании арифметики с плавающей запятой. В этом техническом документе эта тема довольно хорошо освещена.В качестве простой концепции для понимания, если у меня есть числа, представленные как float количества, и я пытаюсь сделать это:

100000000 + 1

результат будет: 100000000 (напишите некоторый коди попробуйте сами)

Это не уникально для графических процессоров, код процессора будет вести себя так же (попробуйте).

Так что для очень больших сокращений мы добираемся до точки (часто) где мы добавляем очень большие числа к гораздо меньшим числам, и результаты не являются точными с точки зрения «чистой математики».

В этом и заключается основная проблема.В вашем коде процессора, когда вы решаете, что правильный результат должен быть 6,1 * n, проблема умножения не ограничивается добавлением больших чисел к маленьким, которые я только что описал, поэтому вы получаете «точный» результат изтот.

Один из способов доказать или обойти это - использовать представление double вместо float.Это на самом деле не полностью устраняет проблему, но подталкивает разрешение к тому, что он может гораздо лучше представить диапазон чисел здесь.

Следующий код в основном имеет это изменение.Вы можете изменить typedef, чтобы сравнить поведение между float и double.

В коде есть несколько других изменений.Ни одна из них не является причиной несоответствия, которое вы засвидетельствовали.

$ cat t18.cu
    #include <iostream>
    #include <math.h>
    #include <stdlib.h>
    #include <iomanip>
    #include <thrust/reduce.h>
    #include <thrust/execution_policy.h>

    #define BLOCK_SIZE 32
    typedef double ft;
    using namespace std;

    __device__ double my_atomicAdd(double* address, double val)
    {
      unsigned long long int* address_as_ull =
                              (unsigned long long int*)address;
      unsigned long long int old = *address_as_ull, assumed;

      do {
        assumed = old;
        old = atomicCAS(address_as_ull, assumed,
                        __double_as_longlong(val +
                               __longlong_as_double(assumed)));

      // Note: uses integer comparison to avoid hang in case of NaN (since NaN != NaN)
      } while (assumed != old);

      return __longlong_as_double(old);
    }
    __device__ float my_atomicAdd(float* addr, float val){
        return atomicAdd(addr, val);
    }

    __global__ void reduce(ft *g_idata, ft *g_odata, int n) {

    __shared__ ft sdata[BLOCK_SIZE];

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

    sdata[threadIdx.x] = (i < n)?g_idata[i]:0;

    __syncthreads();

    for (int s=1; s < blockDim.x; s *=2)
    {
        int index = 2 * s * threadIdx.x;;

        if ((index +s) < blockDim.x)
        {
            sdata[index] += sdata[index + s];
        }
        __syncthreads();
    }


    if (threadIdx.x == 0)
        my_atomicAdd(g_odata,sdata[0]);
    }




    int main(void){

    unsigned int n=pow(10,8);

    ft *g_idata, *g_odata;

    cudaMallocManaged(&g_idata, n*sizeof(ft));
    cudaMallocManaged(&g_odata, sizeof(ft));
    cout << "n = " << n << endl;
    int blockSize = BLOCK_SIZE;
    int numBlocks = (n + blockSize - 1) / blockSize;
    g_odata[0] = 0;
    for(int i=0;i<n;i++){g_idata[i]=6.1;}


    reduce<<<numBlocks, blockSize>>>(g_idata, g_odata, n);
    cudaDeviceSynchronize();


    cout << g_odata[0] << "\t" << (float)n*g_idata[1] << "\t"<< (float)n*g_idata[1]-g_odata[0]<<endl;

    g_odata[0]=thrust::reduce(thrust::device, g_idata, g_idata+n);

    cout << g_odata[0] << "\t" << (float)n*g_idata[1] << "\t"<< (float)n*g_idata[1]-g_odata[0]<<endl;



    cudaFree(g_idata);
    cudaFree(g_odata);

    }
$ nvcc -o t18 t18.cu
$ cuda-memcheck ./t18
========= CUDA-MEMCHECK
n = 100000000
6.1e+08 6.1e+08 0.00527966
6.1e+08 6.1e+08 5.13792e-05
========= ERROR SUMMARY: 0 errors
$