Différents résultats après réduction

Question

J'ai deux algorithmes de réduction, les deux sont de docs.nvidia, donc ils devraient être corrects mais le premier (très efficace) me donne un mauvais résultat. Le deuxième résultat est meilleur mais je m'attendais à une meilleure précision. Y a-t-il une erreur dans les algorithmes ou est-ce que je fais quelque chose de mal?

#include <stdio.h> 
#include <cuda.h> 
#include <stdlib.h> 
#include <math.h> 
#include "cuda_error.h" 

//Lock definition 
#ifndef __LOCK_H__ 
#define __LOCK_H__ 
struct Lock { 
int *mutex; 
Lock(void) { 
CudaSafeCall(cudaMalloc((void**)&mutex, 
sizeof(int))); 
CudaSafeCall(cudaMemset(mutex, 0, sizeof(int))); 
} 
~Lock(void) { 
cudaFree(mutex); 
} 
__device__ void lock(void) { 
while(atomicCAS(mutex, 0, 1) != 0); 
} 
__device__ void unlock(void) { 
atomicExch(mutex, 0); 
} 
}; 
#endif 
//------------------------- 


const int N = 507904; 
const int blockSize = 256; 
int blocks = N/blockSize; 

template <unsigned int blockSize> 
__global__ void reduce(Lock lock, float *g_idata, float *g_odata, int n) 
{ 
     extern __shared__ int sdata[]; 
     unsigned int tid = threadIdx.x; 
     unsigned int i = blockIdx.x*(blockSize*2) + tid; 
     unsigned int gridSize = blockSize*2*gridDim.x; 
     sdata[tid] = 0; 

     while (i < n) 
     { 
      sdata[tid] += g_idata[i] + g_idata[i+blockSize]; 
      i += gridSize; 
     } 

     __syncthreads(); 

     if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); } 
     if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); } 
     if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); } 

     if (tid < 32) 
     { 
      if (blockSize >= 64) sdata[tid] += sdata[tid + 32]; 
      if (blockSize >= 32) sdata[tid] += sdata[tid + 16]; 
      if (blockSize >= 16) sdata[tid] += sdata[tid + 8]; 
      if (blockSize >= 8) sdata[tid] += sdata[tid + 4]; 
      if (blockSize >= 4) sdata[tid] += sdata[tid + 2]; 
      if (blockSize >= 2) sdata[tid] += sdata[tid + 1]; 
     } 

    if (tid == 0) 
    { 
     lock.lock();   
     *g_odata += sdata[0]; 
     lock.unlock(); 
    } 

} 

__global__ void reduction_sum(Lock lock,float *in, float *out, int N) 
{ 
    extern __shared__ float sf_data[]; 
    int tid = threadIdx.x; 
    int i = blockIdx.x * blockDim.x + threadIdx.x; 

    sf_data[tid] = (i<N) ? in[i] : 0; 

    __syncthreads(); 

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

    if (tid == 0) 
    { 
     lock.lock();   
     *out += sf_data[0]; 
     lock.unlock(); 
    } 
} 
//initializing function 
double los() 
{ 
    return (double)rand()/(double)RAND_MAX; 
} 
//------------------------------------------- 


int main (void) 
{ 

//memory allocation 
    float *a; 
    float *dev_a, *dev_blocksum1, *dev_blocksum2; 
    float s1=0, s2=0, spr=0; 

    a = (float*)malloc(N*sizeof(float)); 
    CudaSafeCall(cudaMalloc((void**)&dev_a, N*sizeof(float))); 
    CudaSafeCall(cudaMemset(dev_a, 0, N*sizeof(float))); 
    CudaSafeCall(cudaMalloc((void**)&dev_blocksum1, sizeof(float))); 
    CudaSafeCall(cudaMalloc( (void**)&dev_blocksum2, sizeof(float) ) ); 
    CudaSafeCall(cudaMemset(dev_blocksum1, 0, sizeof(float))); 
    CudaSafeCall(cudaMemset(dev_blocksum2, 0, sizeof(float))); 
//-------------------- 

//drawing, array fill 
    srand(2403); 
    int i; 
    for (i=0; i<N; i++) 
    { 
     a[i]=los(); 
     spr+=a[i]; 
    } 
    printf("CPU sum: %f \n", spr); 
//copy HtoD 
    CudaSafeCall(cudaMemcpy(dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice)); 
//--------------------- 

Lock lock; 

//reduce 
    reduce<blockSize><<<blocks, blockSize, blockSize*sizeof(float)>>>(lock, dev_a, dev_blocksum1, N); 
    CudaSafeCall( cudaMemcpy (&s1, dev_blocksum1, sizeof(float), cudaMemcpyDeviceToHost ) ); 
    printf("GPU sum1: %f \n", s1); 
//----------------------- 

//reduction_sum 
    reduction_sum<<<blocks, blockSize, blockSize*sizeof(float)>>>(lock, dev_a, dev_blocksum2, N); 
    CudaSafeCall( cudaMemcpy (&s2, dev_blocksum2, sizeof(float), cudaMemcpyDeviceToHost ) ); 
    printf("GPU sum2: %f \n", s2); 
//--------------------- 

    return 0; 
} 

$ CPU sum: 253833.515625 
$ GPU sum1: 14021.000000 
$ GPU sum2: 253835.906250 

Answer 1

Il y a un certain nombre de choses à mentionner.

1. Je ne suis pas sûr que votre vérification d'erreur est valide. Vous n'avez pas montré le fichier qui implémente ceci, et quand je cours votre code avec cuda-memcheck, j'obtiens diverses erreurs rapportées. Ils semblent tous être liés à la fonction de verrouillage. Je ne sais pas pourquoi vous utiliseriez la fonction de verrouillage et je ne le recommande pas. Je ne pense pas que ce soit hors de portée quand vous le pensez, en fonction de la façon dont vous l'utilisez. Je suggère d'utiliser atomicAdd à la place, ce qui devrait être plus rapide et plus simple. Au minimum, commentez l'instruction cudaFree() dans votre destructeur.

2. Vous liez une ancienne présentation. Si vous passez en revue un newer version of it alors je pense que vous verrez qu'il recommande maintenant l'utilisation de volatile. Je ne vais pas réécrire votre code entier pour vous, ni résumer ce livre blanc à nouveau, mais si vous ajoutez simplement volatile à la déclaration de la mémoire partagée à des fins de démonstration, il va résoudre le problème résultant.

3. Votre variable de mémoire partagée est déclarée comme int mais vous float quantités résumait. Cela ne fonctionnera pas comme vous le voulez. Vous pouvez le déclarer comme ceci:

   extern __shared__ volatile float sdata[]; 
   

4. Les modifications ci-dessus obtenir le code « travail » pour moi. L'élément restant est l'écart entre les résultats du processeur et du GPU. Je crois que cela est dû au fait que l'ordre des opérations en virgule flottante n'est pas le même sur le CPU (réduction en série) que sur le GPU (réduction parallèle). Étant donné que l'écart se produit dans le 6e chiffre significatif sur une quantité de float, je dirais que cela se situe bien dans la fourchette de vraisemblance pour la comparaison des résultats à virgule flottante. Si vous souhaitez plus de précision, vous pouvez essayer de passer de float à double. En outre, il existe divers documents à virgule flottante que vous souhaitez lire, qui aideront à la compréhension ici, tels que this one et this one.