Podwójna precyzja OpenCL inna niż podwójna precyzja procesora

Question

Programuję w OpenCL, używając karty GeForce GT 610 w systemie Linux. Moje wyniki podwójnej precyzji CPU i GPU nie są spójne. Mogę opublikować część kodu tutaj, ale najpierw chciałbym wiedzieć, czy ktoś inny zmierzył się z tym problemem. Różnica między wynikami podwójnej precyzji GPU i procesora staje się wyraźna, gdy uruchamiam pętle z wieloma iteracjami. W kodzie nie ma nic szczególnego, ale mogę go opublikować tutaj, jeśli ktoś jest zainteresowany. Wielkie dzięki. Oto mój kod. Proszę wymusić __ i złe formatowanie, ponieważ jestem tu nowy :) Jak widać, mam dwie pętle, a mój kod CPU jest zasadniczo prawie identyczną wersją.

#ifdef cl_khr_fp64 
#pragma OPENCL EXTENSION cl_khr_fp64 : enable 
#elif defined(cl_amd_fp64) 
#pragma OPENCL EXTENSION cl_amd_fp64 : enable 
#else 
#error "Double precision floating point not supported by OpenCL implementation." 

#endif

__kernel void simpar(__global double* fp, __global double* fp1, 
    __global double* fp3, __global double* fp5, 
__global double* fp6, __global double* fp7, 
__global double* fp8, __global double* fp8Plus, 
__global double* x, __global double* v, __global double* acc, 
__global double* keBuf, __global double* peBuf, 
unsigned int prntstps, unsigned int nprntstps, double dt 
) { 
unsigned int m,i,j,k,l,t; 
unsigned int chainlngth=100; 
double dxi, twodxi, dxipl1, dximn1, fac, fac1, fac2, fac13, fac23; 
double ke,pe,tke,tpe,te,dx; 
double hdt, hdt2; 
double alpha=0.16; 
double beta=0.7; 
double cmass; 
double peTemp; 
nprntstps=1001; 
dt=0.01; 
prntstps=100; 
double alphaby4=beta/4.0; 
hdt=0.5*dt; 
hdt2=dt*0.5*dt; 
double Xlocal,Vlocal,Acclocal; 
unsigned int global_id=get_global_id(0); 
if (global_id<chainlngth){ 
Xlocal=x[global_id]; 
Vlocal=v[global_id]; 
Acclocal=acc[global_id]; 
for (m=0;m<nprntstps;m++){ 

for(l=0;l<prntstps;l++){ 
       Xlocal =Xlocal+dt *Vlocal+hdt2*Acclocal; 
       x[global_id]=Xlocal; 
       barrier(CLK_LOCAL_MEM_FENCE); 

       Vlocal =Vlocal+ hdt * Acclocal; 
       barrier(CLK_LOCAL_MEM_FENCE); 

      j = global_id - 1; 
      k = global_id + 1; 
      if (j == -1) { 
        dximn1 = 0.0; 
      } else { 
        dximn1 = x[j]; 
      } 
      if (k == chainlngth) { 
        dxipl1 = 0.0; 
      } else { 
        dxipl1 = x[k]; 
      } 
      dxi = Xlocal; 
      twodxi = 2.0 * dxi; 
      fac = dxipl1 + dximn1 - twodxi; 
      fac1 = dxipl1 - dxi; 
      fac2 = dxi - dximn1; 
      fac13 = fac1 * fac1 * fac1; 
      fac23 = fac2 * fac2 * fac2; 
      Acclocal = alpha * fac + beta * (fac13 - fac23); 

      barrier(CLK_GLOBAL_MEM_FENCE); 

      Vlocal += hdt * Acclocal; 
      v[global_id]=Vlocal; 
      acc[global_id]=Acclocal; 
      barrier(CLK_GLOBAL_MEM_FENCE); 
     } 
      barrier(CLK_GLOBAL_MEM_FENCE); 

      tke = tpe = te = dx = 0.0; 
      ke=0.5*Vlocal*Vlocal;//Vlocal*Vlocal; 
      barrier(CLK_GLOBAL_MEM_FENCE); 
      fp6[(m*100)+global_id]=ke; 
      keBuf[global_id]=ke; 
      ke=0.0; 
      barrier(CLK_GLOBAL_MEM_FENCE); 


      j = global_id - 1; 
      k = global_id + 1; 
      if (j == -1) { 
        dximn1 = 0.0; 
      } else { 
        dximn1 = x[j]; 
      } 
      if (k == chainlngth) { 
        dxipl1 = 0.0; 
      } else { 
        dxipl1 = x[k]; 
      } 
      dxi = Xlocal; 
      twodxi = 2.0 * dxi; 
      fac = dxipl1 + dximn1 - twodxi; 
      fac1 = dxipl1 - dxi; 
      fac2 = dxi - dximn1; 
      fac13 = fac1 * fac1 * fac1; 
      fac23 = fac2 * fac2 * fac2; 
      Acclocal = alpha * fac + beta * (fac13 - fac23); 

      barrier(CLK_GLOBAL_MEM_FENCE); 

      Vlocal += hdt * Acclocal; 
      v[global_id]=Vlocal; 
      acc[global_id]=Acclocal; 
      barrier(CLK_GLOBAL_MEM_FENCE); 
     } 
      barrier(CLK_GLOBAL_MEM_FENCE); 

      tke = tpe = te = dx = 0.0; 
      ke=0.5*Vlocal*Vlocal;//Vlocal*Vlocal; 
      barrier(CLK_GLOBAL_MEM_FENCE); 
      fp6[(m*100)+global_id]=ke; 
      keBuf[global_id]=ke; 
      ke=0.0; 
      barrier(CLK_GLOBAL_MEM_FENCE); 
      j = global_id - 1; 
      k = global_id + 1; 
      if (j == -1) { 
        dximn1 = 0.0; 
      } else { 
        dximn1 = x[j]; 
      } 
      if (k == chainlngth) { 
        dxipl1 = 0.0; 
      } else { 
        dxipl1 = x[k]; 
      } 
      dxi = Xlocal; 
      twodxi = 2.0 * dxi; 
      fac = dxipl1 + dximn1 - twodxi; 
      fac1 = dxipl1 - dxi; 
      fac2 = dxi - dximn1; 
      fac13 = fac1 * fac1 * fac1; 
      fac23 = fac2 * fac2 * fac2; 
      Acclocal = alpha * fac + beta * (fac13 - fac23); 

      barrier(CLK_GLOBAL_MEM_FENCE); 

      Vlocal += hdt * Acclocal; 
      v[global_id]=Vlocal; 
      acc[global_id]=Acclocal; 
      barrier(CLK_GLOBAL_MEM_FENCE); 
     } 
      barrier(CLK_GLOBAL_MEM_FENCE); 

      tke = tpe = te = dx = 0.0; 
      ke=0.5*Vlocal*Vlocal;//Vlocal*Vlocal; 
      barrier(CLK_GLOBAL_MEM_FENCE); 
      fp6[(m*100)+global_id]=ke; 
      keBuf[global_id]=ke; 
      ke=0.0; 
      barrier(CLK_GLOBAL_MEM_FENCE); 
    if (global_id ==0){ 
      for(t=0;t<100;t++) 
        tke+=keBuf[t]; 
      } 

      barrier(CLK_GLOBAL_MEM_FENCE); 
      k = global_id-1; 
      if (k == -1) { 
       dx = Xlocal; 
      }else{ 
       dx = Xlocal-x[k]; 
      } 

       fac = dx * dx; 
       peTemp = alpha * 0.5 * fac + alphaby4 * fac * fac; 
       fp8[global_id*m]=peTemp; 
       if (global_id == 0) 
        tpe+=peTemp; 

       barrier(CLK_GLOBAL_MEM_FENCE); 
       cmass=0.0; 
       dx = -x[100-1]; 
       fac = dx*dx; 

       pe=alpha*0.5*fac+alphaby4*fac*fac; 
       if (global_id==0){ 
       fp8Plus[m]=pe; 
       tpe+=peBuf[0]; 
       fp5[m*2]=i; 
       fp5[m*2+1]=cmass; 
       te=tke+tpe; 
       fp[m*2]=m; 
       fp[m*2+1]=te; 

      } 
    barrier(CLK_GLOBAL_MEM_FENCE); 
       //cmass /=100; 
      fp1[(m*chainlngth)+global_id]=Xlocal-cmass; 
      // barrier(CLK_GLOBAL_MEM_FENCE); 
       fp3[(m*chainlngth)+global_id]=Vlocal; 
      // barrier(CLK_GLOBAL_MEM_FENCE); 
      fp7[(m*chainlngth)+global_id]=Acclocal; 

       barrier(CLK_GLOBAL_MEM_FENCE); 
    } 
} 

}

Answer 1

ta jest nieco oczekiwane zachowanie, rzeczywiście.

W starszych procesorach x86 liczby zmiennoprzecinkowe mają długość 80 bitów (Intel "long double"), a obcięte do 64 bitów tylko w razie potrzeby. Gdy jednostki SIMD/instrukcje dla arytmetyki zmiennoprzecinkowej przybyły dla procesorów x86, podwójna precyzja zmiennoprzecinkowa stała się domyślnie 64-bitowa; jednak w zależności od ustawień kompilatora nadal możliwe jest 80-bitowe. Jest dużo do przeczytania na ten temat: Wikipedia: Floating Point.

Sprawdź ustawienia kompilatora dla OpenCL i kodu hosta na zmiennoprzecinkowe "magiczne sztuczki", aby uzyskać lepszą zgodność wyników. Oblicz wartości swoich wartości i sprawdź, czy ten margines błędu jest bezpieczny dla twojej aplikacji.