Двойная точность OpenCL отличается от двойной точности ЦП

Я программирую в OpenCL, используя карту GeForce GT 610 в Linux. Мои результаты двойной точности ЦП и ГП не согласованы. Я могу опубликовать часть кода здесь, но сначала я хотел бы знать, сталкивался ли кто-нибудь еще с этой проблемой. Разница между результатами двойной точности графического процессора и процессора становится заметной, когда я запускаю циклы с большим количеством итераций. На самом деле в коде нет ничего особенного, но я могу опубликовать его здесь, если кому-то интересно. Большое спасибо. Вот мой код. Пожалуйста, извините за __ и плохое форматирование, так как я здесь новичок :) Как видите, у меня есть два цикла, и код моего процессора практически идентичной версии.

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

}


person Newbee    schedule 05.08.2013    source источник
comment
Когда вы говорите «несовместимо», что вы на самом деле имеете в виду? Можете ли вы количественно проиллюстрировать различия (абсолютные и относительные погрешности, количество знаков последнего разряда в согласии)?   -  person talonmies    schedule 05.08.2013
comment
Я выполняю те же (два) вложенных цикла. Внутренний цикл выполняется 100 раз, а затем внешний цикл выводит некоторые результаты. Ниже приведен пример из первой итерации внешнего цикла: Процессор: 0,0000000011 0,0000030832 0,0005244239 0,1400572807 1,5213598941 1,5213598941 0,1400572807 0,0005244239 0,0000030832 0,0000000011 ГПУ: 0,0000000012 0,0000032002 0,0005183133 0,1401775827 1,5249561626 1,5249561626 0,1401775827 0,0005183133 0,0000032002 0,0000000012   -  person Newbee    schedule 05.08.2013
comment
Когда вы говорите о двойной точности ЦП, вы имеете в виду выполнение одного и того же ядра на устройстве ЦП или сравниваете с собственной реализацией? Также убедитесь, что вы не передаете какие-либо флаги компилятору OpenCL, которые могут ослабить строгость IEEE (например, включение быстрых операций mul-add и т. д.). И да, код, безусловно, поможет, чтобы люди могли воспроизводить и проводить свои собственные тесты, без кода лучшее, что мы можем сделать, это строить догадки.   -  person Thomas    schedule 05.08.2013
comment
Извините, я должен был написать ясно. Это нативная однопоточная реализация на ЦП, с которой я сравниваю. Не могли бы вы уточнить, какой это может быть флаг? Это во время постановки ядра в очередь? Спасибо   -  person Newbee    schedule 05.08.2013
comment
@Newbee Это упрощает дело - устройства OpenCL обычно (если не указано иное) работают в строгом соответствии со стандартом IEEE с плавающей запятой. ЦП (в собственном коде) обычно этого не делают, поскольку их часто заставляют либо использовать приближения для повышения производительности, либо работать в 80-битном режиме с плавающей запятой (так называемый расширенный режим), который, как вы можете догадаться, приводит к более точным Результаты. Чтобы быть полностью точным, вам нужно сбросить состояние процессора с плавающей запятой, что немного сложно. Попробуйте вызвать _fpreset() или что-то подобное для вашего языка/платформы перед нативным кодом.   -  person Thomas    schedule 05.08.2013
comment
(на самом деле это большая проблема для многих приложений, особенно сетевого программного обеспечения, которое должно убедиться, что все используемые ЦП действительно выполняют одни и те же вычисления с плавающей запятой, чтобы результаты в конечном итоге не расходились)   -  person Thomas    schedule 05.08.2013
comment
Спасибо большое. Значит, я могу с самого начала вызвать _fpreset() из своего кода на C? Хорошо, я буду читать. Большое спасибо за Вашу помощь.   -  person Newbee    schedule 05.08.2013
comment
Поскольку вы не используете какие-либо трансцендентные функции, должна быть возможность заставить opencl и собственный процессор производить точное совпадение. Если вы можете опубликовать полный пример кода как для opencl, так и для собственного исполнения, я мог бы его отладить. У меня только видеокарта AMD, а не nvidia.   -  person    schedule 05.08.2013
comment
@ScottD, большое спасибо за ваше предложение. Однако, поскольку у нас разные графические процессоры, может быть сложно сопоставить его с графическими процессорами nVidia. С другой стороны, может ли кто-нибудь сказать мне, как использовать _fpreset() из gcc? Я не мог найти никакой информации и не мог определить правильный включаемый файл.   -  person Newbee    schedule 06.08.2013


Ответы (1)


На самом деле это несколько ожидаемое поведение.

На старых процессорах x86 числа с плавающей запятой имеют длину 80 бит (Intel "long double") и усекаются. на 64 бит только при необходимости. Когда модули/инструкции SIMD для арифметики с плавающей запятой появились для процессоров x86, двойная точность с плавающей запятой по умолчанию стала 64-битной; однако 80-битная версия все еще возможна, в зависимости от настроек вашего компилятора. Об этом можно много прочитать: Википедия: Плавающая точка.

Проверьте настройки компилятора для OpenCL и хост-кода на "магические трюки" с плавающей запятой, чтобы добиться лучшего совпадения результатов. Вычислить абсолютную и относительная ошибка ваших значений и проверьте, безопасна ли эта погрешность для вашего приложения.

person Sven    schedule 05.08.2013
comment
Большое спасибо! Я использовал опцию -ffloat-store с gcc, разницы нет. Не могли бы вы предложить какой-либо другой вариант компилятора для gcc, который мог бы соответствовать строгим стандартам IEEE? - person Newbee; 06.08.2013
comment
@Newbee Сводка параметров gcc, связанных с точностью с плавающей запятой: gcc.gnu.org/wiki/FloatingPointMath - person sbabbi; 06.08.2013
comment
Попробуйте использовать -funsafe-math-optimizations -O3 в худшем случае и снизьте до -O0. Кроме того, проверьте параметры переключателя -mfpmath=unit и переключателей -m___ для генерации кода для разных единиц измерения. Я предполагаю, что код OpenCL сильно оптимизирован, поэтому включение небезопасных оптимизаций должно приблизить результаты вашего процессора к результатам OpenCL. - person Sven; 06.08.2013