Differences

This shows you the differences between two versions of the page.

--- etc:common_activities:intel_students_cup:tour2 [2007/11/01 02:47] – kel
+++ etc:common_activities:intel_students_cup:tour2 [2008/01/03 02:32] (current) – external edit 127.0.0.1
@@ Line 1: / Line 1: @@
 ====== Исследование вопроса ======
-пока я выложу свои изыскания. если появятся альтернативы -- сделаем разные странички, или секции -- как будет удобнее.
+  * [[tour2_tasks|Текущие задачи]]
+  * [[algo_details|Детали алгоритма]]
 ===== План =====
@@ Line 153: / Line 155: @@
 ==== Менее очевидные шаги ====
+=== Использование Math Kernel Library ===
+[[http://www.intel.com/cd/software/products/asmo-na/eng/307757.htm|Math Kernel Library]]([[http://www.google.ru/url?sa=t&ct=res&cd=1&url=http%3A%2F%2Fwww.intel.com%2Fsoftware%2Fproducts%2Fmkl%2Ftechtopics%2Fmklman52.pdf&ei=txIpR7fmGIfs0AS7oMj7Aw&usg=AFQjCNGAKhvQ7sI3I9hp7BI6QU5mxa-RUg&sig2=ywkFBVb2LBI5chtUSOhN3Q|Описание]])
+  * библиотека в основном хорошо реализует матричные операции;
+Для исследования Intel Kernel Math Library написал следующее:
-  - ускорить счет синуса. например линейной интерполяцией по таблице:
+Код "глупой" программы, которая вызывает не думая стандартный синус.
+<code cpp>
+#include <stdlib.h>
+#include <math.h>
+int main()
+{
+  double X[32][32];
+  double F[32][32];
+  for(int i=0; i<32; ++i)
+    for(int j=0; j<32; ++j)
+      X[i][j] = rand()%1024;
+  for(int x=0;x<800;++x)
+    for(int y=0;y<600;++y)
+      for(int i=0;i<32;++i)
+        for(int j=0;j<32;++j)
+           F[i][j] = sinf(X[i][j]);
+  return 0;
+}
+</code>
+Код программы, с использованием MKL:
+В нем применил функцию, вычисляющую синус элементов вектора.
+<code cpp>
+#include "mkl.h"
+#include <stdlib.h>
+int main()
+{
+  double X[32][32];
+  double F[32][32];
+  for(int i=0; i<32; ++i)
+    for(int j=0; j<32; ++j)
+      X[i][j] = rand()%1024;
+  for(int x=0;x<800;++x)
+    for(int y=0;y<600;++y)
+      vdSin(32*32,(const double *)X,(double *)F);
+  return 0;
+}
+</code>
+MakeFile:
+<code>
+default: stupid fast
+main.o: main.cpp
+        g++ main.cpp -c -o main.o
+stupid: main.o
+        g++ main.o -o stupid
+clean:
+        rm -f main.o stupid
+fast: imkl_main.o
+        g++ -L/opt/intel/mkl/9.1.023/lib/32 imkl_main.o -lguide -lmkl_p4m -lmkl_ia32 -lm -lirc -o fast
+imkl_main.o: imkl_main.cpp
+        g++ -I/opt/intel/mkl/9.1.023/include -c imkl_main.cpp -o imkl_main.o
+</code>
+Третьим шагом было изменение в "быстрой" программе всех double на float. и вызов функции vsSin
+Результат запуска:
+<code bash>
+make && time ./stupid && time ./fast && time ./floatfast
+real    0m41.288s
+user    0m39.579s
+sys     0m0.143s
+real    0m18.878s
+user    0m18.158s
+sys     0m0.086s
+real    0m7.799s
+user    0m7.412s
+sys     0m0.049s
+</code>
+помимо всего прочего MKL имеет реализацию __одновременного__ вычисления синуса и косинуса в одной функции
+Очевидно, что если в "тупую" программу добавить рядом с вызовом синуса вызов косинуса, то врямя возрастет в два раза, что и произошло при опытной проверке.
+В случае же с MKL, интересней. Далее сравнительные времена выполнения двух программ с применением MKL:
+  - вызываются vsSin и vsCos
+  - вызывается vsSinCos
+<code>
+real    0m16.126s
+user    0m15.261s
+sys     0m0.026s
+real    0m13.289s
+user    0m12.670s
+sys     0m0.028s
+</code>
+=== Использование Intel C Compiler ===
+[[http://www.intel.com/cd/software/products/asmo-na/eng/284132.htm|Компилятор]]([[http://www.intel.com/cd/software/products/asmo-na/eng/compilers/clin/277618.htm|Описание]])
+  * должен позволить автоматически ипользовать SIMD команды (SSE, SSE2...) для оптимизации вычислений в основном в циклах;
+  * возможно подскажет где что можно ещё распараллелить;
+[[icc|Результаты]]
+=== Использование Integrated Performance Primitives ===
+[[http://www.intel.com/cd/software/products/asmo-na/eng/302910.htm|Intel Performance Primitives]]({{intel:sw5r.pdf|Описание}})
+  * как заявляется производительность растёт в том числе и за счёт оптимизации библиотеки под различные модели процессоров;
+  * можно попробовать исопльзовать оттуда не только тригонометрию но и функции работы с изображениями 2D;
+P.S. Менеджер проекта этой библиотеки из Нижнего Новгорода откуда и сам sunset :)
+тестовая программка, а-ля zps:
+<code cpp>
+#include <stdlib.h>
+#ifdef USE_IPP
+#include <ippvm.h>
+#endif
+#include <math.h>
+#include <stdio.h>
+#define COUNT(a) (sizeof(a)/sizeof(*(a)))
+float a[1024];
+float r1[1024];
+float r2[1024];
+int main()
+{
+  for(size_t i=0;i<COUNT(a);++i)
+    a[i]=(drand48()-.5)*20;
+#ifdef USE_IPP
+  for(int i=0;i<800*600;++i)
+    ippsSin_32f_A21(a,r1,COUNT(r1));
+#else
+  for(int i=0;i<800*600;++i)
+    for(size_t j=0;j<COUNT(r1);++j)
+      r2[j]=sinf(a[j]);
+#endif
+#if 0
+  double s=0;
+  for(size_t i=0;i<COUNT(r1);++i)
+  {
+    s+=fabs(r1[i]-r2[i]);
+  }
+  printf("%lg\n",s);
+#endif
+}
+</code>
+<code>
+default: stupid fast
+stupid: main.o
+  g++ main.o -o stupid
+clean:
+  rm -f main.o stupid
+fast: ipp_main.o
+  g++ -L/opt/intel/ipp/5.2/ia32/sharedlib ipp_main.o -lippcore -lippvm -o fast
+ipp_main.o: main.cpp
+  g++ -I/opt/intel/ipp/5.2/ia32/include -DUSE_IPP -c main.cpp -o ipp_main.o
+</code>
+результаты:
+<code>
+$ time ./fast && time ./stupid
+real    0m3.666s
+user    0m3.661s
+sys     0m0.002s
+real    0m37.244s
+user    0m37.095s
+sys     0m0.046s
+</code>
+итого -- ускорение в 10 раз. замена A21 на A11 дает
+<code>
+$ time ./fast && time ./stupid
+real    0m2.975s
+user    0m2.964s
+sys     0m0.007s
+real    0m36.754s
+user    0m36.675s
+sys     0m0.033s
+</code>
+Для тригонометрии в hotspot 2 использовал ippsSin.
 <file>
-diff -ruN src-org/sunset.cpp src/sunset.cpp
+diff -ruN src-org/Makefile src1/Makefile
+--- src-org/Makefile    2007-09-17 17:43:08.000000000 +0400
++++ src1/Makefile       2007-11-02 00:49:53.000000000 +0300
+@@ -56,19 +56,19 @@
+ CC         = gcc
+ CPLAT      =
+ CPROC      =
+-CINC       = -I$(SRC_DIR)
++CINC       = -I$(SRC_DIR) -I/opt/intel/ipp/5.2/ia32/include
+ CDEFS      =
+ COBJ       = -c -o$(OBJ_DIR)/$@
+-CDEFOPT    = -O2
++CDEFOPT    = -g -pg -fprofile-arcs -ftest-coverage
+ COPT       =
+-CFLAGS     =
++CFLAGS     = -O3 -ffast-math -ffinite-math-only -fno-math-errno -funsafe-math-optimizations -fno-trapping-math -march=prescott
+ CFLAGS_ALL = $(CFLAGS) $(CINC) $(CDEFS) $(CDEFOPT) $(CPROC) $(CPLAT)
+ LD         = g++
+ LDPLAT     =
+-LDFLAGS    =
++LDFLAGS    = -g -pg -ax -fprofile-arcs -ftest-coverage -L/opt/intel/ipp/5.2/ia32/sharedlib
+ LDOUTOPT   = -o "$(OUT_DIR)/$(BENCHMARK)"
+-LIBS       = -lm -lc
++LIBS       = -lm -lc -lippcore -lippvm
+ LIBS_ALL   = $(LIBS)
+ endif
+diff -ruN src-org/sunset.cpp src1/sunset.cpp
 --- src-org/sunset.cpp  2007-09-16 12:04:44.000000000 +0400
-+++ src/sunset.cpp  2007-10-24 19:39:04.000000000 +0400
++++ src1/sunset.cpp     2007-11-02 00:34:38.000000000 +0300
-@@ -84,6 +84,27 @@
+@@ -45,6 +45,7 @@
-     return a * (float)k;
+ #include <omp.h>
- }
+ #endif
+ #include "sunset.h"
++#include <ippvm.h>
-+#define SIN_TAB_SZ 1024
+ #define MIN(x,y)    (((x) < (y)) ? (x) : (y))
-+#define PI2 (2*3.141592653589793f)
+ #define MAX(x,y)    (((x) < (y)) ? (y) : (x))
-+float g_sinTab[2*SIN_TAB_SZ+1];
+@@ -747,9 +748,10 @@
-+
-+void fillSinTab()
-+{
-+  for(int i=0;i<2*SIN_TAB_SZ+1;++i)
-+  {
-+    g_sinTab[i]=sinf((i-SIN_TAB_SZ)*PI2/SIN_TAB_SZ);
-+  }
-+}
-+
-+inline float tab_sinf(float v)
-+{
-+  float i=fmodf(v,PI2)*SIN_TAB_SZ/PI2+SIN_TAB_SZ;
-+  size_t idx=i;
-+  float d=i-idx;
-+
-+  return g_sinTab[idx]*(1-d)+g_sinTab[idx+1]*d;
-+}
-+
- /* memory allocation function for two-dimensional arrays */
- void** alloc2D(int M, int N, size_t size)
- {
-@@ -320,6 +341,8 @@
- /*---------------------------------------------------------------------------*/
-     if(iCurFrame == 1)
-     {
-+        fillSinTab();
-+
-         iWaveMeshSize = iWaveHarmNum * iAngleHarmNum;
-         if(iAllocated == 1)
+                 pFlTmp = flArgSin[currentthread].aptr;
-@@ -749,7 +772,7 @@
-                 #pragma ivdep
+-                #pragma ivdep
-                 for(t=0; t<iWaveMeshSize; t++)
+-                for(t=0; t<iWaveMeshSize; t++)
 -                    pFlTmp[t] = (float)sinf(pFlTmp[t]);
-+                    pFlTmp[t] = (float)tab_sinf(pFlTmp[t]);
++                ippsSin_32f_A21(pFlTmp,pFlTmp,iWaveMeshSize);
++                //#pragma ivdep
++                //for(t=0; t<iWaveMeshSize; t++)
++                //    pFlTmp[t] = (float)tab_sinf(pFlTmp[t]);
                  /* initialize the values of derivation */
                  flDerivX = 0.0f;
 </file>
+Результат -- 10.044/кадр, 0.7% отличий. Понижение точности до 11 бит дает 9.139/кадр, 1% отличий. Однако, по-простецки с openmp оно дружить не захотело -- segfault.
+zps.
+Максимальный результат:
+diff -ruN src-org src-ipp
+<code>
+diff src-org/Makefile src-ipp/Makefile
+c59
+< CINC       = -I$(SRC_DIR)
+---
+> CINC       = -I$(SRC_DIR) -I/opt/intel/ipp/5.2/ia32/include
+c62,63
+< CDEFOPT    = -O2
+---
+> #CDEFOPT    = -O2
+> CDEFOPT    = -g
+c65
+< CFLAGS     =
+---
+> CFLAGS     = -O3 -ffast-math -ffinite-math-only -fno-math-errno -funsafe-math-optimizations -fno-trapping-math -march=pentium4m
+c70
+< LDFLAGS    =
+---
+> LDFLAGS    = -g -pg -ax -fprofile-arcs -ftest-coverage -L/opt/intel/ipp/5.2/ia32/sharedlib
+c72
+< LIBS       = -lm -lc
+---
+> LIBS       = -lm -lc -lippcore -lippvm -lguide -lipps -lippm
+diff src-org/sunset.cpp src-ipp/sunset.cpp
+c48,51
+<
+---
+> #include <ippvm.h>
+> #include <ipps.h>
+> #include <ippm.h>
+> #include <iostream>
+c310,311
+<                 float   OT, KX1, KY1;
+---
+>                 float   *OT, *KX1, *KY1;
+>                 float   kx1, ky1, ot; // zps
+a344,346
+>             free(OT);
+>             free(KX1);
+>             free(KY1);
+c364,366
+<
+---
+>         OT           = (float*)malloc(iWaveHarmNum * sizeof(float));
+>         KX1          = (float*)malloc(iWaveHarmNum * sizeof(float));
+>         KY1          = (float*)malloc(iWaveHarmNum * sizeof(float));
+a372
+>
+a668,669
+>     ippsMulC_32f(flOmega, flTime, OT, iWaveHarmNum);
+>     float xxx[iAngleHarmNum];
+c675
+<         private(currentthread, OT, KX1, KY1) \
+---
+>         private(currentthread, OT, KX1, KY1, kx1, ky1, ot) \
+a692
+>
+c745,749
+< */
+---
+> *///              std::cerr << "NKMAX: "<< NKMAX << std::endl;
+>
+>                 ippsMulC_32f(flK, flDecartX[i][j], KX1, NKMAX);
+>                 ippsMulC_32f(flK, flDecartY[i][j], KY1, NKMAX);
+>
+,738c752,755
+<                     OT  = flOmega[t] * flTime;
+<                     KX1 = flK[t] * flDecartX[i][j];
+<                     KY1 = flK[t] * flDecartY[i][j];
+<
+---
+>                     kx1 = KX1[t];
+>                     ky1 = KY1[t];
+>                     ot =  OT[t];
+>                     int len = t * iAngleHarmNum;
+,745c758,761
+<                         iSinIndex1 = t * iAngleHarmNum + l;
+<                         flArgSin[currentthread].aptr[iSinIndex1] = OT -
+<                             KX1 * flAzimuthCosFi[l] - KY1 * flAzimuthSinFi[l] +
+<                             flRandomPhase[t*iAngleHarmNum + l];
+<                     } /* end for l */
+---
+>                         iSinIndex1 = len + l;
+>                         flArgSin[currentthread].aptr[iSinIndex1] = ot -
+>                            kx1*flAzimuthCosFi[l] - ky1*flAzimuthSinFi[l] + flRandomPhase[iSinIndex1];
+>                    } /* end for l */
+c763,776
+<
+---
+> #if 0
+>                 float * dest;
+>                 for(t = 0; t < NKMAX; t++)
+>                 {
+>                     dest = &flArgSin[currentthread].aptr[t*iAngleHarmNum];
+>                     kx1 = KX1[t];
+>                     ky1 = KY1[t];
+>                     ippmLComb_vv_32f(flAzimuthCosFi, 4,kx1, flAzimuthSinFi, 4,ky1, xxx, 4, iAngleHarmNum);
+>                     ippsSubCRev_32f(xxx,OT[t],dest, iAngleHarmNum);
+>
+>                 } /* end for t */
+>                 dest = NULL;
+> #endif
+> //                ippsAdd_32f_I(flRandomPhase, flArgSin[currentthread].aptr, iWaveMeshSize);
+,752c779,783
+<                 #pragma ivdep
+<                 for(t=0; t<iWaveMeshSize; t++)
+<                     pFlTmp[t] = (float)sinf(pFlTmp[t]);
+---
+> 		ippsSin_32f_A21(pFlTmp,pFlTmp,iWaveMeshSize);
+>
+>                 //#pragma ivdep
+>                 //for(t=0; t<iWaveMeshSize; t++)
+>                 //    pFlTmp[t] = (float)sinf(pFlTmp[t]);
+a790,801
+> #if 0
+>                 {
+>                   const float *pr[]={flAmplitudeX,flAmplitudeY};
+>                   float __r[2];
+>
+>                   ippmDotProduct_vav_32f_L(pr,0,4,pFlTmp,4,__r,iWaveMeshSize,2);
+>
+>                   flDerivX = __r[0];
+>                   flDerivY = __r[1];
+>                 }
+> #endif
+> #if 1
+a811
+> #endif
+Только в src-ipp: .sunset.cpp.swp
+</code>
+результат:
+<code bash>
+Frame 1 of 16 ... frame time 8.534
+Frame 2 of 16 ... frame time 8.508
+Frame 3 of 16 ... frame time 8.167
+Frame 4 of 16 ... frame time 8.538
+Frame 5 of 16 ... frame time 8.467
+Frame 6 of 16 ... frame time 8.302
+Frame 7 of 16 ... frame time 8.565
+Frame 8 of 16 ... frame time 8.425
+Frame 9 of 16 ... frame time 8.549
+Frame 10 of 16 ... frame time 8.398
+Frame 11 of 16 ... frame time 8.425
+Frame 12 of 16 ... frame time 8.595
+Frame 13 of 16 ... frame time 8.419
+Frame 14 of 16 ... frame time 8.312
+Frame 15 of 16 ... frame time 8.471
+Frame 16 of 16 ... frame time 8.417
+=================================
+Timing:
+        Total time is 135.098 sec., average frame time is 8.444 sec.
+Correctness check:
+        Max RGB difference is 11.
+        Number of different color pixels is 3564 (0.7%).
+</code>
-Даже такая простецкая реализация дает время 23.724/кадр, или ускорение в 1.5 раза. При этом количество отличающихся пикселов растет по сравнению с оригинальным алгоритмом незначительно. Если интерполяцию заменить просто поиском ближайшего узла таблицы, точность падает значительно (условие 5% отличных пикселов не выполняется).
+{{:intel:result_kel.txt|Результаты работы на конфигурации P4(2,6)Windows}}:\\
-Полагаю, более тонкая (ассемблерная?) реализация синуса даст еще небольшой прирост скорости.
+Лучший из полученных результатов:
+<code bash>
+Frame 1 of 16 ... frame time 3.591
+Frame 2 of 16 ... frame time 3.492
+Frame 3 of 16 ... frame time 3.497
+Frame 4 of 16 ... frame time 3.499
+Frame 5 of 16 ... frame time 3.485
+Frame 6 of 16 ... frame time 3.492
+Frame 7 of 16 ... frame time 3.494
+Frame 8 of 16 ... frame time 3.492
+Frame 9 of 16 ... frame time 3.523
+Frame 10 of 16 ... frame time 3.491
+Frame 11 of 16 ... frame time 3.496
+Frame 12 of 16 ... frame time 3.495
+Frame 13 of 16 ... frame time 3.494
+Frame 14 of 16 ... frame time 3.492
+Frame 15 of 16 ... frame time 3.503
+Frame 16 of 16 ... frame time 3.495
+=================================
+Timing:
+        Total time is 56.036 sec., average frame time is 3.502 sec.
-Откомпилировать другими компиляторами. если времена для одного и того же кода будут заметно отличаться -- смотреть в ассемблер. или просто использовать другой компилятор и дальше (.
+Correctness check:
+        Max RGB difference is 20.
+        Number of different color pixels is 5437 (1.1%).
+</code>
+Достигнут:\\
+  указанием ключа fast при компилляции
+  заменой дублирующих вычислений в циклах в очевидных местах
+  выносом вычисления синусов и косинусов (ipp) полярных координат
-В части испльзования инструментов intel можно применить следующее:\\
+==== Совсем неочевидные/исследовательские шаги ====
-- Провести компилляцию с их [[http://www.intel.com/cd/software/products/asmo-na/eng/284132.htm|компиллятором]]([[http://www.intel.com/cd/software/products/asmo-na/eng/compilers/clin/277618.htm|Описание]]):\\
+  * [[tab_sin|Табличный синус]] (плавающие числа) -- тупиковая ветвь.
-  а) должен позволить автоматически ипользовать SIMD команды (SSE, SSE2...) для оптимизации вычислений в основном в циклах\\
+  * Целочисленная реализация алгоритма расчета вектора нормали -- выигрыша по скорости нет, по точности -- проигрыш. Тупиковая ветвь.
-  б) возможно подскажет где что можно ещё распараллелить\\
+  * [[harm_reduction|Уменьшение количества гармоник]] -- тупиковая ветвь.
-- Для тригонометрии использовать [[http://www.intel.com/cd/software/products/asmo-na/eng/302910.htm|Intel Performance Primitives]]({{intel:sw5r.pdf|Описание}}):\\
-  а) как заявляется производительность растёт в том числе и за счёт оптимизации библиотеки под различные модели процессоров\\
-  б) можно попробовать исопльзовать оттуда не только тригонометрию но и функции работы с изображениями 2D\\
-P.S. Менеджер проекта этой библиотеки из Нижнего Новгорода откуда и сам sunset :)\\
-- Попытаться использовать [[http://www.intel.com/cd/software/products/asmo-na/eng/307757.htm|Math Kernel Library]]([[http://www.google.ru/url?sa=t&ct=res&cd=1&url=http%3A%2F%2Fwww.intel.com%2Fsoftware%2Fproducts%2Fmkl%2Ftechtopics%2Fmklman52.pdf&ei=txIpR7fmGIfs0AS7oMj7Aw&usg=AFQjCNGAKhvQ7sI3I9hp7BI6QU5mxa-RUg&sig2=ywkFBVb2LBI5chtUSOhN3Q|Описание]]):\\
-  а) Как более глубокий вариант оптимизации с учётом того что библиотека восновном хорошо реализует матричные операции