cuda dot produit de deux vecteurs ne fonctionnant pas pour N> = 369

Question

J'ai mis en œuvre un produit vectoriel dot comme ci-dessous. Il est compilé avec CUDA 7.5 avec compute_20,sm_20 et const int THREADS_PER_BLOCK = 16;.

Cela se produit pour les flotteurs et les doubles.

Cela fonctionne jusqu'à n=368, mais au-delà, le résultat est incorrect. Je me demandais si le problème venait de mon code d'implémentation ou des valeurs que j'utilise (voir le deuxième code, les initialisations), par ex. peut être l'ajout au-delà de n=368 introduit des erreurs à virgule flottante (cela peut être étrange car la même erreur se produit à la fois pour float et doubles).

int divUp(int total, int grain) { return (total+grain-1)/grain; } 

__device__ __forceinline__ double 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))); 
    } 
    while(assumed!=old); 
    return __longlong_as_double(old); 
} 

__device__ __forceinline__ float atomicAdd(float* address, float val) 
{ 
    unsigned int *ptr = (unsigned int *)address; 
    unsigned int old, newint, ret = *ptr; 
    do { 
     old = ret; 
     newint = __float_as_int(__int_as_float(old)+val); 
    } while((ret = atomicCAS(ptr, old, newint)) != old); 

    return __int_as_float(ret); 
} 

template<typename T> 
__global__ void vecdotk(const T* a, const T* b, const int n, T* c) 
{ 
    __shared__ T temp[THREADS_PER_BLOCK]; 
    int x = threadIdx.x+blockIdx.x*blockDim.x; 
    if(x==0) c[0] = 0.0; 
    if(x<n) {temp[threadIdx.x] = a[x]*b[x]; 
    } 
    else temp[threadIdx.x] = 0.0; 
    __syncthreads(); 

    if(0==threadIdx.x) 
    { 
     T sum = 0.0; 
     for(int j=0; j<THREADS_PER_BLOCK; ++j) 
     { 
      sum += temp[j]; 
     } 
     atomicAdd(c, sum); 
    } 
} 

template<typename T> 
void dot(const T* a, const T* b, const int n, T* c) 
{ 
    dim3 block(THREADS_PER_BLOCK); 
    dim3 grid(divUp(n, block.x), 1); 
    vecdotk<T><<<grid, block>>>(a, b, n, c); 
    cudaSafeCall(cudaGetLastError()); 
}; 

J'utilise les deux vecteurs hôtes suivants pour remplir les matrices de dispositif d'entrée (qui a je ne suis pas représentés car ils font partie d'une plus grande bibliothèque). Fondamentalement, je veux calculer la somme de la série carrée c.-à-

// fill host vectors a and b 
for(int i=0; i<n; ++i) 
{ 
    h_vec_a[i] = i+1;//__mat_rand(); 
    h_vec_b[i] = i+1;//__mat_rand(); 
} 

Answer 1

Cela ne fonctionnera pas:

if(x==0) c[0] = 0.0; 

Il n'y a aucune garantie (en CUDA) que le fil 0 se dirige d'abord ou que cette ligne sera exécutée avant que les autres threads arrivent à un point quelconque du code. Vous devrez initialiser c[0] avant de lancer ce noyau. Dans le cas contraire, certains threads peuvent faire leur ajout atomique à c, puis, un peu plus tard, le thread 0 peut initialiser c [0] à zéro.

En outre, CUDA fournit déjà une version float de atomicAdd, il n'y a aucune raison pour que vous fournissiez la vôtre. Et, l'exécution de threads de 16 threads ne vous donnera pas de bonnes performances (je recommande d'utiliser simplement la fonction produit CUBLAS.) Avec le correctif pour c[0] (supprimer cette ligne de code, et initialiser c[0] avant le noyau) votre code fonctionne correctement pour moi:

$ cat t372.cu 
#include <stdio.h> 

const int n = 2048; 
#ifdef USE_DOUBLE 
typedef double mytype; 
#else 
typedef float mytype; 
#endif 
const int THREADS_PER_BLOCK = 16; 

int divUp(int total, int grain) { return (total+grain-1)/grain; } 
#if 0 
__device__ __forceinline__ double 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))); 
    } 
    while(assumed!=old); 
    return __longlong_as_double(old); 
} 

__device__ __forceinline__ float atomicAdd(float* address, float val) 
{ 
    unsigned int *ptr = (unsigned int *)address; 
    unsigned int old, newint, ret = *ptr; 
    do { 
     old = ret; 
     newint = __float_as_int(__int_as_float(old)+val); 
    } while((ret = atomicCAS(ptr, old, newint)) != old); 

    return __int_as_float(ret); 
} 
#endif 
template<typename T> 
__global__ void vecdotk(const T* a, const T* b, const int n, T* c) 
{ 
    __shared__ T temp[THREADS_PER_BLOCK]; 
    int x = threadIdx.x+blockIdx.x*blockDim.x; 
    //if(x==0) c[0] = 0.0; 
    if(x<n) {temp[threadIdx.x] = a[x]*b[x]; 
    } 
    else temp[threadIdx.x] = 0.0; 
    __syncthreads(); 

    if(0==threadIdx.x) 
    { 
     T sum = 0.0; 
     for(int j=0; j<THREADS_PER_BLOCK; ++j) 
     { 
      sum += temp[j]; 
     } 
     atomicAdd(c, sum); 
    } 
} 

template<typename T> 
cudaError_t dot(const T* a, const T* b, const int n, T* c) 
{ 
    dim3 block(THREADS_PER_BLOCK); 
    dim3 grid(divUp(n, block.x), 1); 
    vecdotk<T><<<grid, block>>>(a, b, n, c); 
    cudaDeviceSynchronize(); 
    return cudaGetLastError(); 
}; 

int main(){ 

    mytype *h_vec_a, *h_vec_b, *d_vec_a, *d_vec_b, *h_c, *d_c; 
    int bs = n*sizeof(mytype); 
    h_vec_a = (mytype *)malloc(bs); 
    h_vec_b = (mytype *)malloc(bs); 
    h_c = (mytype *)malloc(sizeof(mytype)); 
    cudaMalloc(&d_vec_b, bs); 
    cudaMalloc(&d_vec_a, bs); 
    cudaMalloc(&d_c, sizeof(mytype)); 
// fill host vectors a and b 
    for(int i=0; i<n; ++i) 
    { 
    h_vec_a[i] = i+1;//__mat_rand(); 
    h_vec_b[i] = i+1;//__mat_rand(); 
    } 
    h_c[0] = 0; 
    cudaMemcpy(d_vec_a, h_vec_a, bs, cudaMemcpyHostToDevice); 
    cudaMemcpy(d_vec_b, h_vec_b, bs, cudaMemcpyHostToDevice); 
    cudaMemcpy(d_c, h_c, sizeof(mytype), cudaMemcpyHostToDevice); 
    dot(d_vec_a, d_vec_b, n, d_c); 
    cudaMemcpy(h_c, d_c, sizeof(mytype), cudaMemcpyDeviceToHost); 
    mytype test_val = 0; 
    for (int i=0; i < n; i++) 
    test_val += h_vec_a[i] * h_vec_b[i]; 
    printf("GPU result: %f, CPU result: %f\n", h_c[0], test_val); 

} 
$ nvcc -arch=sm_20 -o t372 t372.cu 
nvcc warning : The 'compute_20', 'sm_20', and 'sm_21' architectures are deprecated, and may be removed in a future release (Use -Wno-deprecated-gpu-targets to suppress warning). 
$ cuda-memcheck ./t372 
========= CUDA-MEMCHECK 
GPU result: 2865411584.000000, CPU result: 2865411072.000000 
========= ERROR SUMMARY: 0 errors 
$ 

la différence numérique dans les 3 derniers chiffres est due aux limites de float, non dues à une erreur dans le code. Si vous changez votre initialisation pour initialiser chaque vecteur à tous les 1, par exemple, vous obtiendrez un résultat exact dans ce cas.

Encore une fois, pour des raisons de performances, un certain nombre de critiques pourraient être formulées à l'encontre de votre code. Si vous voulez un produit scalaire rapide, je suggère d'utiliser CUBLAS.