Dot Produkt mit einem CUDA Kernel für große Vektorgrößen liefert falsche Ergebnisse

Question

Ich versuche einen CUDA Kernel zu implementieren, der das Skalarprodukt zweier Vektoren berechnet. Für kleine Vektorgrößen funktioniert der Code richtig und ich bekomme die richtigen Ergebnisse, aber für größere Vektoren wird nur falsch berechnet. Ich war die Umsetzung drei verschiedene Möglichkeiten, um das Punktprodukt zu berechnen:

• serielle Version (gerade in C++ ohne Optimierung)
• CUDA Kernel
• CUBLAS Version

Mein Haupt in der cpp-Datei sieht wie folgt aus:

float *h_x,*h_y; 
float res1=0.0, res2=0.0, res3=0.0; 

h_x=(float*)malloc(Main::N*sizeof(float)); random_ints_Vec(h_x); 
h_y=(float*)malloc(Main::N*sizeof(float)); random_ints_Vec(h_y); 

double serialTimer; 
double cublasTimer; 
double cudaTimer; 

res1=serial_dotProd(h_x,h_y,&serialTimer); 
res2=cublas_dotProd(h_x,h_y,&cublasTimer); 
res3=cuda_dotProd(h_x,h_y,&cudaTimer);  

free(h_x); free(h_y); 

Serien Version:

float Main::serial_dotProd(float* x, float* y, double* time){ 
std::clock_t start; 
start = std::clock(); 

float res1=0.0; 
for (int i=0;i<Main::N;++i) { 
    res1+=x[i]*y[i]; 
} 

*time= (std::clock() - start)/(double) CLOCKS_PER_SEC; 
return res1;} 

CUDA Version:

float Main::cuda_dotProd(float *h_x,float *h_y,double* time){ 
float *d_x,*d_y,*d_res, *h_res; 
h_res=(float*)malloc(Main::BLOCKS_PER_GRID*sizeof(float)); 

size_t bfree, afree, total; 
cudaMemGetInfo(&bfree,&total); 

cudaMalloc((void**) &d_x,Main::N*sizeof(float)); 
cudaMalloc((void**) &d_y,Main::N*sizeof(float)); 
cudaMalloc((void**) &d_res,Main::BLOCKS_PER_GRID*sizeof(float)); 
cudaCheckErrors("cuda malloc fail"); 

cudaMemGetInfo(&afree,&total); 
std::cout<<" > memory used for cuda-version:"<< (bfree -afree)/1048576<< "MB ("<<total/1048576 <<"MB)" <<"\n"; 

cudaMemcpy(d_x,h_x,Main::N*sizeof(float),cudaMemcpyHostToDevice); 
cudaMemcpy(d_y,h_y,Main::N*sizeof(float),cudaMemcpyHostToDevice); 

std::clock_t start; 
start = std::clock(); 
DotProdWrapper(d_x,d_y,d_res,(Main::N+Main::THREADS_PER_BLOCK-1)/Main::THREADS_PER_BLOCK,Main::THREADS_PER_BLOCK,Main::N); 

*time= (std::clock() - start)/(double) CLOCKS_PER_SEC; 

cudaMemcpy(h_res,d_res,Main::BLOCKS_PER_GRID*sizeof(float),cudaMemcpyDeviceToHost); 

float c= 0; 
for (int i=0;i<Main::BLOCKS_PER_GRID;++i){ 
    c+=h_res[i]; 
} 
cudaFree(d_x); 
cudaFree(d_y); 
cudaFree(d_res);  

free(h_res); 
return c;} 

CUDA Kernel:

__global__ void DotProd(float* x, float* y, float* scalar,unsigned long int N){ 
    extern __shared__ float cache[]; 

    int tid = threadIdx.x+ blockIdx.x*blockDim.x; 
    int cacheIndex = threadIdx.x; 

    float temp=0; 
    while (tid<N){ 
     temp+=x[tid]*y[tid]; 
     tid +=blockDim.x*gridDim.x; 
    } 
    cache[cacheIndex]=temp; 
    __syncthreads(); 

    int i=blockDim.x/2; 
    while(i!=0){ 
     if (cacheIndex<i) 
      cache[cacheIndex]+=cache[cacheIndex+i]; 
     __syncthreads(); 
     i/=2; 
    } 
    if(cacheIndex==0) 
     scalar[blockIdx.x]=cache[cacheIndex]; 
} 

CUBLAS Version:

float Main::cublas_dotProd(float *h_x,float *h_y, double* time){ 
    float *d_x,*d_y; 
    float *res; 
    float result=0.0; 
    cublasHandle_t h; 
    cublasCreate(&h); 
    cublasSetPointerMode(h, CUBLAS_POINTER_MODE_DEVICE); 

    size_t bfree, afree, total; 
    cudaMemGetInfo(&bfree,&total); 

    cudaMalloc((void**) &d_x,Main::N*sizeof(float)); 
    cudaMalloc((void**) &d_y,Main::N*sizeof(float)); 
    cudaMalloc((void **)(&res), sizeof(float)); 
    cudaCheckErrors("cublas malloc fail"); 

    cudaMemGetInfo(&afree,&total); 

    cudaMemcpy(d_x, h_x, Main::N*sizeof(float), cudaMemcpyHostToDevice); 
    cudaMemcpy(d_y, h_y, Main::N*sizeof(float), cudaMemcpyHostToDevice); 

    cublasSetVector(Main::N,sizeof(float),h_x,1,d_x,1); 
    cublasSetVector(Main::N,sizeof(float),h_y,1,d_y,1); 


    std::clock_t start; 
    start = std::clock(); 

    cublasSdot(h,Main::N,d_x,1,d_y,1,res); 

    *time= (std::clock() - start)/(double) CLOCKS_PER_SEC; 

    cudaMemcpy(&result, res, sizeof(float), cudaMemcpyDeviceToHost); 

    cudaFree(d_x); 
    cudaFree(d_y); 
    cudaFree(res); 
    return result; 
} 

Die Ergebnisse, die ich mit den verschiedenen Einstellungen nach der Berechnung erhalten werden aufgelistet:

• N = 204.800, THREADS_PER_BLOCK: 512, BLOCKS_PER_GRID: 400 serial_dotProd = 4.15369e + 06; cublas_dotProd = 4.15369e + 06; cuda_dotProd = 4.15369e + 06
• N = 20.480.000, THREADS_PER_BLOCK: 512, BLOCKS_PER_GRID: 40000 serial_dotProd = 4.04149e + 08; cublas_dotProd = 4.14834e + 08; cuda_dotProd = 4.14833e + 08

Ich weiß nicht, warum, aber nach einer gewissen Größe meiner Vektoren bekomme ich nur falsches Ergebnis. Die Vektoren passen in das SDRAM und der gemeinsame Speicher für jeden Block hat auch genug Platz, um den Speicher zuzuordnen. Vielen Dank im Voraus.

Answer 1

Diese Frage ist so oft gestellt worden, dass Nvidia an entire section of the CUDA Floating Point and IEEE 754 Anleitung dazu gewidmet hat. Es wird auch kurz in der CUDA C Best Practices Guide erwähnt.

Die kurze Erklärung dafür ist zweifach:

• Im Gegensatz zu den entsprechenden exakten mathematischen Operationen, arithmetische Gleitkomma-Operationen wegen der Rundungsfehler beteiligt nicht assoziativ sind. Dies bedeutet, dass das Umordnen der Summierung von einer geraden seriellen Summe in eine Baumstruktur, die für die parallele Ausführung geeignet ist, das Ergebnis geringfügig ändert (mehr mit zunehmender Anzahl von summierten Werten). Zufälligerweise ergibt die Baumanordnung in den meisten Fällen auch ein Ergebnis, das näher an der exakten mathematischen Summe liegt als die sequentielle Summe.

• Der CUDA-Compiler neigt dazu, bei der Verwendung von Fused Multiply-Add (FMA, eine Multiply-Though-Add-Operation, bei der die Zwischenrundungsstufe weggelassen wird) aggressiver zu sein. Wiederum tendiert das mathematisch korrekte Ergebnis dazu, näher an dem unter Verwendung von FMA erhaltenen Ergebnis zu liegen.

So ist die wahrscheinliche Antwort ist, dass die mit CUDA erhaltenen Ergebnisse sind wahrscheinlich näher an das genaue Ergebnis als eine einfache serielle CPU-Version (weshalb ich Sie gebeten, das Experiment durchzuführen wieder erhöhte Präzision verwendet wird).