Оптимизация умножения вектора 1хnC на матрицу nRxnC в SSE2

Тема в разделе "WASM.BEGINNERS", создана пользователем Skevalt, 15 апр 2007.

  1. Skevalt

    Skevalt New Member

    Публикаций:
    0
    Регистрация:
    15 апр 2007
    Сообщения:
    3
    Доброго времени суток всем! В ходе работы над нейронными сетями встал вопрос об оптимизации алгоритма обучения однойслойного персептрона. В частности, наибольшее время занимают две функции - это умножение матрицы nR x nC на вектор nC и сложение двух матриц такой же размерности.
    Важно использовать тип данных double(ошибка копиться за 100-200 эпох существенная). Необходимо реализовать функции на SSE2. LjДо этого все было в классах и алгоритм школьный- так что заморочек не было;)
    Попытался написать алгоритм умножения, который пришел в голову на Си:
    Код (Text):
    1. const int nR = 10;
    2. const int nC = 25;
    3.  
    4.  
    5. _MM_ALIGN16 double  *inVector,
    6.             *outVector,
    7.             *lentMatrix;
    8. //---------------------------------------------------------------------------
    9. void InitializiVectorAndMatrix(int nRow, int nCol, double *ma, double *vec)
    10. {
    11.    int i, j;
    12.    for(i = 0; i < nRow; i++)
    13.       for(j = 0; j < nCol; j++)
    14.       {
    15.          ma[i*nCol + j] = (double) (rand()%1000)/1000.0;
    16.          vec[j] = (double) (rand()%1000)/1000.0;
    17.       }
    18. }
    19. //---------------------------------------------------------------------------
    20. void MatrixMulVectorAdvanced(int nRow, int nCol, double *ma, double *vecIn, double *vecOut)
    21. {
    22.    int i, j;
    23.  
    24.    double x11, x12, x21, x22, x31, x32, x41, x42, y11, y12, y21, y22, y31, y32;
    25.  
    26.    int intNumCol_2 = 2*(nCol/2);
    27.    int intNumCol_6 = 6*(nCol/6);
    28.  
    29.    int shift = 0;
    30.    int nshift = 0;
    31.  
    32.    i = 0;
    33.    do{
    34.       x41 = 0.0;
    35.       x42 = 0.0;
    36.  
    37.       j = 0;
    38.      shift = i*nCol;
    39.       while(j < intNumCol_6){
    40.          nshift = shift + j;
    41.          x11 = ma[nshift + 0];
    42.          x12 = ma[nshift + 1];
    43.  
    44.          x21 = ma[nshift + 2];
    45.          x22 = ma[nshift + 3];
    46.  
    47.          x31 = ma[nshift + 4];
    48.          x32 = ma[nshift + 5];
    49.  
    50.          y11 = vecIn[j + 0];
    51.          y12 = vecIn[j + 1];
    52.  
    53.          y21 = vecIn[j + 2];
    54.          y22 = vecIn[j + 3];
    55.  
    56.          y31 = vecIn[j + 4];
    57.          y32 = vecIn[j + 5];
    58.  
    59.          x41 += x11*y11;
    60.          x42 += x12*y12;
    61.  
    62.          x41 += x21*y21;
    63.          x42 += x22*y22;
    64.  
    65.          x41 += x31*y31;
    66.          x42 += x32*y32;
    67.  
    68.          j+= 6;
    69.       }
    70.  
    71.       while(j < intNumCol_2){
    72.          nshift = shift + j;
    73.          x11 = ma[nshift + 0];
    74.          x12 = ma[nshift + 1];
    75.  
    76.          y11 = vecIn[j + 0];
    77.          y12 = vecIn[j + 1];
    78.  
    79.          x41 += x11*y11;
    80.          x42 += x12*y12;
    81.            j += 2;
    82.       }
    83.  
    84.       x41 += x42;
    85.  
    86.       while(j < nCol){
    87.          x11 = ma[shift + j];
    88.          y11 = vecIn[j];
    89.  
    90.          x41 += x11*y11;
    91.            j++;
    92.       }
    93.  
    94.       vecOut[i] = x41;
    95.  
    96.    }while(++i < nRow);
    97. }
    А потом что-то написал на assemblere...Куча ошибок и совершенно неоптимизарованно, но только час занимаюсь ассемблером. Не могли бы вы мне помочь разобраться с этим кодом?
    Код (Text):
    1. void MatrixMulVectorAdvancedSSE2(int nRow, int nCol, double *ma, double *vecIn, double *vecOut)
    2. {
    3.     int intNumCol_2 = 2*(nCol/2);
    4.     int intNumCol_6 = 6*(nCol/6);
    5.     __asm{
    6.         pushad
    7.  
    8.         xor eax, eax            ;//Counter for rows in matrix
    9.         xor ebx, ebx            ;//Counter for cols in matrix and rows in vector
    10.         mov esi, ma             ;//Load pointer to source matrix
    11.         mov edi, vecIn          ;//Load pointer to source vector   
    12.        
    13.         Loop_By_Rows:           ;//Cycle for rows of matrix
    14.             cmp eax, nRow
    15.             jnl End_Loop_By_Rows
    16.             xorpd xmm6, xmm6
    17.             xor ebx, ebx
    18.  
    19.             mov edx, intNumCol_6
    20.             Loop_By_Col_Step_6:
    21.                 cmp ebx, edx
    22.                 jnl Loop_By_Col_Step_2
    23.                 mov  ecx, nCol
    24.                 imul ecx, eax
    25.                 add  ecx, ebx  
    26.                 movupd xmm0, [esi + ecx + 0h]
    27.                 movupd xmm1, [esi + ecx + 10h]
    28.                 movupd xmm2, [esi + ecx + 20h]
    29.  
    30.                 movupd xmm3, [edi + ecx + 0]
    31.                 movupd xmm4, [edi + ecx + 10h]
    32.                 movupd xmm5, [edi + ecx + 20h]
    33.  
    34.                 mulpd  xmm0, xmm3
    35.                 mulpd  xmm1, xmm4
    36.                 mulpd  xmm2, xmm5
    37.  
    38.                 addpd  xmm6, xmm0
    39.                 addpd  xmm6, xmm1
    40.                 addpd  xmm6, xmm2
    41.  
    42.                 add ebx, 6  
    43.             jmp Loop_By_Col_Step_6
    44.  
    45.             Loop_By_Col_Step_2:
    46.                 mov edx, intNumCol_2
    47.                 cmp ebx, edx
    48.                 jnl End_Loop_By_Col_Step_2
    49.                 mov  ecx, nCol
    50.                 imul ecx, eax
    51.                 add  ecx, ebx
    52.                 movupd xmm0, [esi + ecx + 0]
    53.                 movupd xmm3, [edi + ecx + 0]
    54.  
    55.                 mulpd  xmm0, xmm3
    56.                 addpd  xmm6, xmm0
    57.                 add ebx, 2
    58.             jmp Loop_By_Col_Step_2
    59.  
    60.             End_Loop_By_Col_Step_2:
    61.             movupd   xmm1, xmm6
    62.             shufpd   xmm1, xmm1, 3h
    63.             addsd    xmm6, xmm1
    64.  
    65.             Loop_By_Col_Step_1:
    66.                 mov edx, nCol
    67.                 cmp ebx, edx
    68.                 jnl End_Loop_By_Col_Step_1
    69.                 mov  ecx, nCol
    70.                 imul ecx, eax
    71.                 add  ecx, ebx
    72.                 movsd xmm0, [esi + ecx + 0]
    73.                 movsd xmm3, [edi + ecx + 0]
    74.  
    75.                 mulpd  xmm0, xmm3
    76.                 addpd  xmm6, xmm0
    77.  
    78.                 inc ebx
    79.             jmp Loop_By_Col_Step_1
    80.  
    81.             End_Loop_By_Col_Step_1:
    82.             mov edx, vecOut
    83.             mov ecx, eax
    84.             imul ecx, 8
    85.             movsd [edx + ecx], xmm6
    86.        
    87.         inc eax
    88.         jmp  Loop_By_Rows
    89.  
    90.         End_Loop_By_Rows:
    91.         popad
    92.     }
    93. }
     
  2. TermoSINteZ

    TermoSINteZ Синоби даоса Команда форума

    Публикаций:
    2
    Регистрация:
    11 июн 2004
    Сообщения:
    3.578
    Адрес:
    Russia
    Skevalt
    Для начала глянь, что генерирует Си компиллер интела со вкл оптимизацией по скорости.
     
  3. Skevalt

    Skevalt New Member

    Публикаций:
    0
    Регистрация:
    15 апр 2007
    Сообщения:
    3
    У меня траффика не хватает скачать Win C++ компилер:) Дизассемблировал то, что получается в MSVS компиляторе - но там и не пахнет SSE2. С другой стороны, осталось немного - просто передвижение по массиву не работает в асмвском коде корректно - что-то я упускаю из вида.
    Так как все равно писать dll-йку для Buildera, то надо доразбираться с кодом так...Уже крыша едет от ассемблера;)
     
  4. Stiver

    Stiver Партизан дзена

    Публикаций:
    0
    Регистрация:
    18 дек 2004
    Сообщения:
    812
    Адрес:
    Germany
    Skevalt

    Можно в принципе скачать Intel Math Kernel Library и посмотреть, как это там сделано.
     
  5. leo

    leo Active Member

    Публикаций:
    0
    Регистрация:
    4 авг 2004
    Сообщения:
    2.542
    Адрес:
    Russia
    Skevalt
    Во-во ;) Разворот основного цикла на 6 даже при правильной реализации является избыточным, а в приведенном сишном варианте выглядит просто глупо - спрашивается куда долго и упорно грузятся 12 переменных, если регистров всего 8 ?!!
    Суть разворота подобных циклов заключается не столько в уменьшении накладных расходов на управление циклом (т.к.целочисленные операци могут выполняться параллельно с fpu), сколько в скрытии латентности зависимых операций fadd\addpd за счет параллельного накопления нескольких частичных сумм. Поэтому вариант SSE с точки зрения оптимизации по скорости тоже реализован не правильно, т.к. он копит только одну сумму, а нужно копить как минимум две, как в сишном варианте. И потом использование невыравненной загрузки movupd вместо выравненной movapd может свести на нет все преимущества SSE перед fpu
    Вариант с накоплением 2-х сумм
    Код (Text):
    1.   push esi
    2.   push edi
    3.   push ebx
    4.   mov esi,ma
    5.   mov ebx,vecOut
    6.   mov eax,nCol
    7.   mov edx,eax
    8.   and edx,3
    9.   xor eax,edx
    10.   shl eax,3
    11.   neg eax
    12.   mov esp,nRow  ;!!! используется esp, в Си можно заменить на локальную переменную
    13. Loop_By_Rows:
    14.   mov edi,vecIn
    15.   sub edi,eax
    16.   mov ecx,eax
    17.   xorps xmm0,xmm0
    18.   xorps xmm1,xmm1
    19. @@:
    20.   movapd xmm2,[esi]        ;на P4 movupd до 50% хуже movapd
    21.   movapd xmm3,[esi+16]
    22.   mulpd xmm2,[edi+ecx]
    23.   mulpd xmm3,[edi+ecx+16]
    24.   addpd xmm0,xmm2
    25.   addpd xmm1,xmm3
    26.   add esi,32
    27.   add ecx,32
    28.   jnz @B
    29.  
    30.   test edx,2
    31.   jz @F
    32.   movapd xmm2,[esi]
    33.   mulpd xmm2,[edi]
    34.   add esi,16
    35.   add edi,16
    36.   addpd xmm0,xmm2
    37.  
    38. @@:
    39.   test edx,1
    40.   jz @F
    41.   movsd xmm3,[esi]
    42.   mulsd xmm3,[edi]
    43.   add esi,8
    44.   addsd xmm1,xmm3
    45. @@:
    46.   addpd xmm0,xmm1
    47.   movhlps xmm1,xmm0
    48.   addsd xmm1,xmm0
    49.   movq [ebx],xmm1      ;или movsd, movlpd - без разницы
    50.   add ebx,8
    51.   sub esp,1                 ;!!!
    52.   jnz Loop_By_Rows
    53.  
    54.   lea esp,[ebp-4*3]      ;!!!
    55.   pop ebx
    56.   pop edi
    57.   pop esi
    58.   ;ret
     
  6. Skevalt

    Skevalt New Member

    Публикаций:
    0
    Регистрация:
    15 апр 2007
    Сообщения:
    3
    Тут немного другое, когда я писал на С, то подразумевал x11 и х12 страшую и младшую части регистра xmm, таким образом, я гружу всего в три регистра 6 чисел из матрицы, а в другие три регистра 6 чисел из вектора - вот откуда загрузка "12" переменных:) Получается, что за раз сначала перемножается 6 пар чисел, затем, в следующем цикле перемножаются по 2 пары чисел, и потом уже идет работа с одиночными парами.
    Развернул так, чтоб поэффективней использовать сразу все регистры(на мой взгляд) - ведь минимум происходит параллельное сложение двух чисел в старшей и младшей частях регистра.
    А в сишном тоже одна сумма копиться:)
    А это я откровенно ступил - долго искал, как бы выравнивать данные в MSVS (в Builder просто делал указатель с адресом кратным 16) чтобы использовать MOVAPD и в конечном итоге забыл поменять в самом коде
    Спасибо большое за отзывы и код - буду разбираться:)
     
  7. leo

    leo Active Member

    Публикаций:
    0
    Регистрация:
    4 авг 2004
    Сообщения:
    2.542
    Адрес:
    Russia
    Угу, но компилятор не понял твоих замыслов и затолкал все переменные в стек :)))

    Кроме регистров нужно еще эффективно использовать порты запуска и исполнительные блоки. При длинной серии мувов блок load\AGU пыхтит от натуги, а вычислительные блоки отдыхают, хотя при чередовании операций они могли бы работать параллельно с загрузкой

    В твоей интерпретации может и одна, но фактически x41 и x42 это две суммы ;)

    Только учти, что при нечетном nCol половина строк все равно не будет выравнена на 16 и использование movdqa приведет к ошибке - нужны доп.извращения на случай нечетного nCol
     
  8. asmfan

    asmfan New Member

    Публикаций:
    0
    Регистрация:
    10 июл 2006
    Сообщения:
    1.004
    Адрес:
    Abaddon
  9. Ustus

    Ustus New Member

    Публикаций:
    0
    Регистрация:
    8 авг 2005
    Сообщения:
    834
    Адрес:
    Харьков
    asmfan
    Которое до сих пор ни на одном проце эффективно не реализовано :dntknw:
    Разве что подождать аэмдэшной барселоны - они там чем-то подобным пугали... :):):)