OpenCL 双精度与 CPU 双精度不同

我正在 Linux 中使用 GeForce GT 610 卡进行 OpenCL 编程。我的CPU和GPU双精度结果不一致。我可以在这里发布部分代码，但我首先想知道是否有其他人遇到过这个问题。当我运行多次迭代的循环时，GPU 和 CPU 双精度结果之间的差异变得很明显。该代码确实没有什么特别的，但如果有人感兴趣，我可以将其发布在这里。多谢。这是我的代码。请原谅 __ 和错误的格式，因为我是新来的:) 如您所见，我有两个循环，我的 CPU 代码本质上几乎是相同的版本。

#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);
  }
 }

}

实际上，这在某种程度上是预期的行为。

在较旧的 x86 CPU 上，浮点数的长度为 80 位（Intel 的“长双” http://en.wikipedia.org/wiki/Long_double），并且仅在需要时截断为 64 位。 当浮点运算的 SIMD 单元/指令到达 x86 CPU 时，浮点双精度默认变为 64 位；但是，80 位仍然是可能的，具体取决于您的编译器设置。关于这一点，有很多值得阅读的内容：维基百科：浮点 http://en.wikipedia.org/wiki/Floating_point.

检查 OpenCL 的编译器设置and浮点“魔术”上的主机代码，以获得更好的结果一致性。计算absolute http://mathworld.wolfram.com/AbsoluteError.html and 相对误差 http://mathworld.wolfram.com/RelativeError.html您的值并检查此误差范围对于您的应用程序是否安全。