xmm регистры и SSE2

[vitaly333]

Есть 8 чисел двойной точности: A1,A2,B1,B2,C1,C2,C3,C4;
Нужно вычислить (используя 128 битовые xmm регистры и минимальным кол-вом IO-операций между кэшем и регистрами):
С1 = C1 + A1*B1;
C2 = C2 + A1*B2;
C3 = C3 + A2*B1;
C4 = C4 + A2*B2;

Если в xmm0 записать {A1,A1} а в xmm1 - {B1,B2}
Как выполнить операцию xmm3 = xmm0*xmm1; Регистр xmm3 используется для временного хранения произведений A*B.
Меня интересует несколько вопросов:
1) Как будет быстрее: записывать сразу {A1,A1} и в xmm0 и в xmm3 регистры. Т.е. в кэш придется лазить 4 раза. Или загрузить {A1,A1} в xmm0, а потом скопировать данные оттуда в xmm3 регистр и выполнить умножение xmm3 = xmm3*xmm1;
Если последний вариант быстрее, то как выполнить копирование регистр-регистр(не обращаясь к кэшу), используя SSE2 инструкцию.
2) Как правильно и быстро загрузить 2 double числа в реистр? 1-ое число загрузить в младшие 64 разряда xmm-регистра, 2 число - в старшие 64 разряда этого же региста. Или как то упаковать 2 double числа в одно 128-битное или копировать их в double[2] массив и потом загружать его одной инструкцией?

Сообщение отредактировано: vitaly333 - 03.10.09, 11:58

[Pavia]

vitaly333
Думать надо.

C1,C2,C3,C4 - надо хранить в регистрах 2 регистра

Чтение скорее всего лучше\надо сразу чтобы A1,A2,B1,B2, читались соответственно в 2 регистра. А дальше уже раскидывать по регистрам и множить.

[vitaly333]

Цитата

Чтение скорее всего лучше\надо сразу чтобы A1,A2,B1,B2, читались соответственно в 2 регистра. А дальше уже раскидывать по регистрам и множить.

Я немного не корректно написал (хотел упростить задачу). На самом деле есть не 8 а 14 чисел: A1,A2,B1,B2,B3,B4,C1,C2,C3,C4,C5,C6,C7,C8.
И нужно посчитать:
С1 = C1 + A1*B1;
C2 = C2 + A1*B2;
C3 = C3 + A1*B3;
C4 = C4 + A1*B4;

С5 = С5 + A2*B1;
С6 = С6 + A2*B2;
С7 = С7 + A2*B3;
С8 = С8 + A2*B4;

Это можно представть себе как умножение вектора-столбца на вектор-строку:

|C1 C2 C3 C4| = |C1 C2 C3 C4| +  |A1|  X  |B1 B2 B3 B4|
|C5 C6 C7 C8|   |C5 C6 C7 C8|    |A2|

Под x86 у нас в распоряжении есть 8 векторных регистров xmm0-xmm7 так ведь. 4 регистра (xmm4-xmm7) уходят под хранение С элементов. 1 (xmm3) для временного хранения произведений.Остается 3 (xmm0-xmm2). Допустим (как вы говорите) загружаем в xmmo = (A1,A2), в xmm1 = (B1,B2), в
xmm2 = (B3,B4);


Алгоритм :
1. Вычислить С1,С6:
1.1 Копируем xmm0->xmm3
1.2 Вычисляем (A1*B1,A2*B2): xmm3 = xmm3*xmm1;
1.3 Сохраняем результат в xmm4(C1,C6): xmm4 = xmm4 + xmm3
2. Вычислить С3,С8:
2.1 Копируем xmm0->xmm3
2.2 Вычисляем (A1*B3,A2*B4) xmm3 = xmm3*xmm2;
2.3 Сохраняем результат в xmm5(C3,C8): xmm5 = xmm5 + xmm3
3. Вычислить C5,C2:
3.1 Копируем xmm0->xmm3, так чтобы в xmm3 было (A2,A1)
3.2 Вычисляем (A2*B1,A1*B2): xmm3 = xmm3*xmm1;
3.3 Сохраняем результат в xmm6(C5,C2): xmm6 = xmm6 + xmm3
4. Вычислить С7,С4:
4.1 Копируем xmm0->xmm3, так чтобы в xmm3 было (A2,A1)
4.2 Вычисляем (A2*B3,A1*B4): xmm3 = xmm3*xmm2;
4.3 Сохраняем результат в xmm7(C7,C4): xmm7 = xmm7 + xmm3

Это вы имелли ввиду?
Меня интересует как выполнить копирование регистр->регистр минуя кэш (шаги 1.1,2.1 и особенно 3.1 и 4.1) и будут ли это копирование быстрее, чем загружать в xmm3 A1,A2 -коэффициенты прямо из кэша?
И как максимально быстро загружать сразу два double в регистр(например xmmo = (A1,A2)) ?
И вообще является ли такой алгоритм оптимальным?

Добавлено 03.10.09, 13:37

Цитата

Мы можем что-то сказать об их расположении или выравнивании? Или это твоя прграммы, и следовательном мы можем как угодно задать их положение прямо в исходнике? Или об их положении ничего вообще нельзя сказать?

Хранятся вот так :

|C1 C2 C3 C4| = |C1 C2 C3 C4| +  |A1|  X  |B1 B2 B3 B4|
|C5 C6 C7 C8|   |C5 C6 C7 C8|    |A2|

Но, это не значит A1,A2 можно загрузить используя адрес массива A. (A1,A2) это лишь его часть. На самом деле A,B,C - это матрицы N*N размерномти.
А приведенный выше алгоритм это один его самый внутренний цикл. Сами массивы A,B,C могу выровнять по границе 16 байт.

Сообщение отредактировано: vitaly333 - 03.10.09, 20:42

[leo]

Цитата vitaly333 @ 03.10.09, 13:29

Я немного не корректно написал (хотел упростить задачу).
Это можно представть себе как умножение вектора-столбца на вектор-строку

Если ты о задачке умножения больших матриц, то так и следует говорить. В этом случае огромный выигрыш дает транспонирование матрицы B - в первую очередь за счет улучшения условий кэширования, во вторую за счет более эффективного использования SSE2, т.к. в этом случае умножается строка A на строку B. Для повышения эффективности можно за один цикл просчитывать 4 элемента C, производя умножение одной строки A на 4 строки B - тут и регистров с головой хватает и все дубово-линейно без всяких анпаклов\шафлов.

  ;eax = *A; edx = *B1; ebx = *B2; edi = *B3; esi = *B4
  ;ecx = N*8; eax=eax+ecx; edx=edx+ecx; ... esi=esi+ecx
  ;neg ecx
  ;xmm0 = xmm1 = xmm2 = xmm3 = 0 - накапливаемые суммы для C1..C4
loop1:
  movapd xmm4,[eax+ecx]  ;A1
  movapd xmm5,xmm4
  mulpd  xmm4,[edx+ecx]  ;B1
  movapd xmm6,xmm5
  mulpd  xmm5,[ebx+ecx]  ;B2
  movapd xmm7,xmm6
  mulpd  xmm6,[edi+ecx]  ;B3
  mulpd  xmm7,[esi+ecx]  ;B4
  addpd xmm0,xmm4
  addpd xmm1,xmm5
  addpd xmm2,xmm6
  addpd xmm3,xmm7
  add ecx,16
  jnz loop1
;здесь нужно сложить верхние и нижние половинки xmm0..xmm3
;и прибавить их к С1..С4

[vitaly333]

Цитата

Как я мог иметь это ввиду, если ты изначально сформулировал условие другой задачи?

Извините, но "Это вы имелли ввиду?" - это я не вам написал. Это адресовалось Pavia

Цитата

Я бы задумался сначала об его правильности. Например, как ты выбираешь пары вообще не понятно

Я выбираю пары так, чтобы минимизировать загрузки из кэша. Алгоритм, который я выше привел требует на один цикл 3(6) загрузок из кэша на 16 арифм. операций (8 умножений и 8 сложений) + 4 операции копирования из регистра в регистр. При другом порядке вычислений элементов матрицы С их может быть больше.
На самом деле не важно в каком порядке (сначала С1, С2, С3...и.т.д. или сначала С1, С6 потом С3, С8....и.т.д.) вычисляются элементы матрицы С (с точки зрения правильности и корректности вычислений), а важно это для уменьшения запросов в кэш.

Цитата

Зачем 4? Зачем их вообще хранить на регистрах? Если уж экономить регистры, то стоит вспомнить, что инструкция addpd может работать с памятью, т.е. вообще регистры для C не нужны, только для временных переменных при вычислении A*B.

Эти 4 регистра нужны для накапливания произведений A*B, которые потом выгружаются в коэффициенты матрицы С. A регистры не нужно экономить, а нужно задействовать их все. Чем больше регистров - тем выше производительность алгоритма. Нужно лишь организовать вычисления так чтобы меньше лазить в кэш-память.

Цитата

Если ты о задачке умножения больших матриц, то так и следует говорить. В этом случае огромный выигрыш дает транспонирование матрицы B - в первую очередь за счет улучшения условий кэширования, во вторую за счет более эффективного использования SSE2, т.к. в этом случае умножается строка A на строку B. Для повышения эффективности можно за один цикл просчитывать 4 элемента C, производя умножение одной строки A на 4 строки B - тут и регистров с головой хватает и все дубово-линейно без всяких анпаклов\шафлов.

Да именно о них. Я написал вам недавно письмо В ПМ как раз по этой теме, но вы почему то так мне и не ответили. Возможно не дошло. Поэтому я и решил написать здесь. Я не стал писать изначально про умножения матриц потому что блочный алгоритм и так сложный..а тут ещё и регистры и SSE2. Хотел максимально упростить задачу, рассмотрев только внутренний цикл. Но вижу что не совсем у меня это получилось. Придется всё таки описать полностью задачу. Вот часть из моего письма к вам (немного измененная):
Я сейчас занимаюсь анализом блочных алгоритмов умножения матриц и вот написал тут один такой алгоритм:

// Эта функция выполняет блочное умножение матриц вида C = C+A*B
// где,N - размерность матриц
// A,B,C - вещественные квадратные матрицы двойной точности, размерностью [N*N]
// bsize - размер  блока, который должен быть кратен N
double** blockMult(double**A, double** B,double**C, int N,int bsize){
 
        int n,m,iter,i,j,k;
        int ibeg,iend,jbeg,jend,kbeg,kend;
 
        double r,r0,r1,r2,r3; // для эмуляции регистров, накапливающих суммы произведений
        double A0=0.0,B0=0.0,B1=0.0; //
        transpose(B,N); // транспонирование матрицы B
                // исходные матрицы A,B,C мысленно разбиваются на блоки, размером [bsize x bsize]
                // образуя решетку,составленную из этих блоков
        int gridSize = N/bsize; // размер решетки блоков
                // n-ый и m-ный циклы позволяют обойти каждый блок решетки  
                // цикл по блокам строк
        for (n=0; n<gridSize; n++){
            ibeg = n*bsize;
            iend = ibeg+bsize;
                  // цикл по блокам столбцов
            for (m=0; m<gridSize; m++){
                jbeg = m*bsize;
                jend = jbeg+bsize;
                //prefetch(C,ibeg,jbeg,bsize,8); // загрузка блока С[n,m] в кэш
                for (iter=0; iter<gridSize; iter++) {
                    kbeg = iter*bsize;
                    kend = kbeg+bsize;
                    //prefetch(A,ibeg,kbeg,bsize,8); // загрузка блока A[n,iter] в кэш
                    //prefetch(B,jbeg,kbeg,bsize,8); // загрузка блока B[m,iter] в кэш
                    // след. 3 цикла выполняют перемножение блоков:
                     //С[n,m] =  C[n,m]+A[n,iter]*B[m,iter]
                      for (i=ibeg; i<iend; i+=2){
                        for (j=jbeg; j<jend; j+=2){
                            r0 = 0.0;
                            r1 = 0.0;
                            r2 = 0.0;
                            r3 = 0.0;
                            for (k=kbeg; k<kend;k++){
                             // на 2-арифм.операции - 2 загрузки кэш-регистры
                            //C[i][j] += A[i][k]*B[j][k];
                               // алгоритм ниже на
                               //  8 ариф. операций требует
                                // 4 загрузки из кэша в регистры
                                A0 = A[i][k];
                                B0 = B[j][k];
                                B1 = B[j+1][k];
                                r = A0*B0;
                                r0 +=r;
                                r = A0*B1;
                                r1+=r;
                                A0 = A[i+1][k];
                                r = A0*B0;
                                r2 +=r;
                                r = A0*B1;
                                r3+=r;
 
                            }
                            C[i][j] +=r0;
                            C[i][j+1]+=r1;
                            C[i+1][j]+=r2;
                            C[i+1][j+1]+=r3;
                        }
                    }
                }
            }
 
        }
 
 
        return C;
 
 
    }
 
 // Эта функция позволяет загрузить блок матрицы в кэш, прочитав первый double
 // из каждой кэш-линейки
 // ipos,jpos - указывают на позицию блока(верхний левый угол) в матрице A
 // linesize - размер кэш-линии/размер элемента
 void inline prefetch(double** A,int ipos,int jpos,int bsize,int linesize){
 
       int i,j;
       double temp;
       int iend = ipos + bsize;
       int jend = jpos + bsize;
       for (i=ipos;i<iend;i++)
           for (j=jpos;j<jend;j+=linesize)
            temp = A[i][j];
 
   }
 
// Эта функция осуществляет транспонирование матрицы
 void transpose(double** A, int n){
        double tmp = 0.0;
        for (int i=0;i<n;i++){
            for (int j=0;j<i;j++){
                tmp = A[i][j];
                A[i][j] = A[j][i];
                A[j][i] = tmp;
            }
        }
    }

Алгоритм работает правильно и использует блокирование кэша и регистров. Работает довольно быстро (около 50% от пик. производительности, не считая паралельность ядер), но это ещё не всё что можно выжать из него.
Странно, но почему то не работает префетч, который я написал. Пользы от него совершенно никакой. Что с ним, что без него результат один и тот же. Для проверки брал (N = 2016 bsize = 48 linesize = 8).
Размер блока bsize подбирал так, чтобы 3 блока, размером bsize x bsize умещались в кэш L1 и bsize был кратен N и linesize. Именно в подобранных выше цифрах выполняется это условие. Для справки, размер кэша L1 на моем процессоре = 64KB и размер кэш-линейки 64 байта (8 double).

2-ая "странность" заключается в том, что нет никакого эффекта от SSE2 при включенной соответствующей опции в копиляторе. Проверял на GNU G++ 4.4 , MSVS С++ 8 и INTEL C++ v 11.1.Результат один и тот же от SSE2 никакого ускорения нет, а хотя должно быть, особенно на Intel. Может это из-за того что у меня AMD Athlon? Хотя он поддерживает SSE2 и SSE3. При компиляции интелом указывал разные ключи SSE2-оптимизации, в том числе и \QaxW, но результата никакого.

Поэтому я решил написать свою собственную реализацию SSE2, но я к сожалению не знаком с Ассемблером.
Вот накидал тут примерный алгоритм как всё должно выглядеть:

double** blockMultSSE2(double**A, double** B,double**C, int N,int bsize){
 
     int n,m,iter,i,j,k;
     int ibeg,iend,jbeg,jend,kbeg,kend;
 
 
     // для хранения 2-х коэффициентов матрицы A
     double xmm00,xmm01;
     // для хранния 4-х коээфициентов матрицы B
     double xmm10,xmm11,xmm20,xmm21;
     // переменные для эмуляции векторных SIMD-регистров
     // для накапливания произведений
     double xmm40,xmm41,xmm50,xmm51,xmm60,xmm61,xmm70,xmm71;
     // переменные в которые сохранется промежуточный результат произведения A*B
     double xmm30,xmm31;
    
 
 
 
     transpose(B,N);
 
     int gridSize = N/bsize;
        for (n=0; n<gridSize; n++){
            ibeg = n*bsize;
            iend = ibeg+bsize;
            for (m=0; m<gridSize; m++){
               jbeg = m*bsize;
               jend = jbeg+bsize;
                //prefetch(C,ibeg,jbeg,bsize,8); // загрузить блок С[n,m] в кэш
                for (iter=0; iter<gridSize; iter++) {
                  kbeg = iter*bsize;
                  kend = kbeg+bsize;
                      // загрузить  блок A[n,iter] в кэш
                     //prefetch(A,ibeg,kbeg,bsize,8);
                      // загрузить блок B[m,iter] в кэш
                      //prefetch(B,jbeg,kbeg,bsize,8);
                     // след. 3 цикла выполняют перемножение блоков:
                     //С[n,m] =  C[n,m]+A[n,iter]*B[m,iter]
                     for (i=ibeg; i<iend; i+=2){
                       for (j=jbeg; j<jend; j+=4){
                        // обнулить регистры для накапливания произведений xmm4-xmm7
                         xmm40=0.0;  // C1
                         xmm41=0.0;  // C6
                         xmm50=0.0;  // C3
                         xmm51=0.0;  // C8
                         xmm60=0.0;  // C5
                         xmm61=0.0;  // C2
                         xmm70=0.0;  // C7
                         xmm71=0.0;  // C4
                         // далее внутренний цикл  
                        for (k=kbeg; k<kend;k++){
                            // C[i][j] += A[i][k]*B[j][k]; // обычный алгоритм строка-строка
                         // Ниже я реализую алгоритм который я привел в 4 посте
                         // Он позволяет за один цикл просчитывать 8 элементов (частично конечно потому что в блоках)
                         // Сперва нужно загрузить все используемые коэффициенты матриц A и B из кэша
                         // Загрузить коэффициенты матрицы A из кэша в регистр xmm0:
                             xmm00 = A[i][k];   // A1
                             xmm01 = A[i+1][k]; // A2
                         // Загрузить коэффициенты матрицы B из кэша в регистры xmm1,xmm2:
                             xmm10 = B[j][k];   // B1
                             xmm11 = B[j+1][k]; // B2
                             xmm20 = B[j+2][k]; // B3
                             xmm21 = B[j+3][k]; // B4
                         // итого: 3 (или 6 смотря как загружать) операций перемещения кэш-регистры
                         // Теперь сам алгоритм из поста №4: опускаем шаги копирования регистр-регистр, т.к.
                         // здесь используются обычные переменные, а не регистры.
                        // 1. Вычислить С1,С6                                                
                        // 1.2 Перемножить xmm0*xmm1 и результат записать в xmm3 для временного хранения
                             // два след. выражения должны ввыполнятся за один такт    
                             xmm30 = xmm00*xmm10; // A1*B1
                             xmm31 = xmm01*xmm11; // A2*B2
                        // 1.3 Результат добавить к регистру xmm4
                              // так же за один такт
                            xmm40 +=xmm30; // С1 = С1 + A1*B1
                            xmm41 +=xmm31; // C6 = C6 + A2*B2
                        // 2. Вычислить С3,С8
                        // 2.2 Перемножить xmm0*xmm2 и результат поместить в xmm3 для временного хранения
                            xmm30 = xmm00*xmm20; // A1*B3
                            xmm31 = xmm01*xmm21; // A2*B4
                        // 2.3 Результат добавить к регистру xmm5
                            xmm50 +=xmm30; // C3 = C3 + A1*B3;
                            xmm51 +=xmm31; // C8 = C8 + A2*B4;
                        // 3. Вычислить С5,С2
                        // 3.2 Перемножить xmm0*xmm1 и результат поместить в xmm3 для временного хранения
                            xmm30 = xmm01*xmm10; // A2*B1
                            xmm31 = xmm00*xmm11; // A1*B2
                        // 3.3 Результат добавить к регистру xmm6
                            xmm60 +=xmm30;     // C5 = C5 + A2*B1
                            xmm61 +=xmm31;     // C2 = C2 + A1*B2
                        // 4. Вычислить С7,С4
                        // 4.2 Перемножить xmm0*xmm2 и результат поместить в xmm3 для временного хранения
                            xmm30 = xmm01*xmm20; // A2*B3
                            xmm31 = xmm00*xmm21; // A1*B4
                        // 4.3 Результат добавить к регистру xmm7
                            xmm70 +=xmm30; // C7 = C7 + A2*B3
                            xmm71 +=xmm31; // C4 = C4 + A1*B4
                                }                        
                        // выгрузить регистры в переменные матрицы С
                                C[i][j]  +=xmm40; // C1
                                C[i][j+1]+=xmm61; // C2
                                C[i][j+2]+=xmm50; // C3
                                C[i][j+3]+=xmm71; // C4
                                C[i+1][j]+=xmm60; // C5
                                C[i+1][j+1]+=xmm41; //C6
                                C[i+1][j+2]+=xmm70; //C7
                                C[i+1][j+3]+=xmm51; //C8
                            }
                        }
                    }
                }
 
            }
   return C;
 }

Вот такой вот алгоритм.Не тривиальный конечно. Думаю с комментариями будет понятен.Вы не могли бы помочь написать внутренний цикл и prefetch на инлайн-ассемблере? Очень нужно.

Скорее всего он будет требовать 6 загрузок из кэша/на один внутренний цикл, т.к. элементы блоков A и B читаются из столбцов, т.е. разных участков памяти, а это значит что 2 double нужно будет загружать 2-мя инструкциями (MOVSD и MOVHPD) - это 6 тактов (если не ошибаюсь) + 16 операций (8 умножений + 8 сложений) - это 8 тактов + 4 копирования регистр-регистр это 4 такта.
Итого 18 тактов процессора (если всё правильно посчитал) для того чтобы "частично" (в пределах блока) вычислить 8 элементов.
Впринципе можно добиться того чтобы на один внутренний цикл было 3 загрузки из кэша. Есть два пути:
1) Создать рабочие массивы double A[2] double B[4] и перед тем как начинать вычисления копированить в них нужные элементы A и B а потом адреса этих двух массивов загружать в регистры.А из них уже, используя MOVAPD, загружать нужные коэффициенты.
2) Сделать так чтобы внутри каждого блока доступ к элементам осуществлялся строго последовательно, т.е. грубо говоря блок это есть линейный массив размера [bsize*bsize], но тогда A** уже не подойдет и нужна спец. структура данных типа BlockMatrix, которая содержит массив указателей на блоки.
Как лучше будет сделать?
И какие ещё оптимизации можно использовать (не считая многопоточности), чтобы достичь пиковой производительности этого алгоритма?

Сообщение отредактировано: vitaly333 - 04.10.09, 09:34

[leo]

Цитата vitaly333 @ 03.10.09, 20:36

Да именно о них. Я написал вам недавно письмо В ПМ как раз по этой теме, но вы почему то так мне и не ответили. Возможно не дошло

Дошло. Только, во-первых, это тема не для ПМ, во-вторых, ты валишь в одну кучу множество вопросов, да еще предлагаешь свои решения\реализации, в которых с ходу "без поллитры" не разберешься
Поэтому нужно, во-первых, четко\кратко сформулировать основную задачу - умножение двух больших матриц A и B с предварительным транспонированием матрицы B. "Большая" матрица означает, что она не умещается в кэш L2.
Во-вторых, разделить два основных вопроса:
1) Оптимизация внутреннего цикла - распараллеливание вычислений как за счет векторизации (SSE), так и за счет сокращения числа загрузок из памяти на 1 цикл (использование одного вх.значения для вычисления нескольких выходных). По сути - это и есть вопрос данной темы.
2) Имеет ли смысл с точки зрения условий кэширования использовать блочную обработку по сравнению с перемножением в лоб, какие должны быть размеры блоков С и A,B (далеко не факт, что нужно юзать одинаковые и к тому же квадратные блоки для всех матриц). Это задачка скорее чисто алгоритмическая, т.к. ее можно формализовать не прибегая к знанием тонкостей работы процессора. Чтобы не валить все в одну кучу, видимо лучше выделить эту задачку в отдельную тему или в этом разделе или в Алгоритмах, ну или по крайней мере вернуться к ней в этой же теме, но после завершения "разбирательств" с первой задачей

[vitaly333]

leoЭто конечно. Я поэтому и хотел упростить задачу, сформулировав её изначально по другому, чтобы не получалось каши.

Цитата

Имеет ли смысл с точки зрения условий кэширования использовать блочную обработку по сравнению с перемножением в лоб

Имеет.Если вы под "перемножением в лоб" подрузамеваете алгортим строка-строка, то он работает эффективно только до тех пор пока 3 строки умещаются к кэш. Блочный же алгортим эффективен на любой размерности матриц A,B,C.

Цитата

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

Видимо да. Как получу ответ на 1 вопрос можно будет перенести тему в Алгоритмы и там обсудить второй.

[leo]

Цитата vitaly333 @ 04.10.09, 09:30

Я поэтому и хотел упростить задачу, сформулировав её изначально по другому, чтобы не получалось каши.

Угу, только в итоге ты ее сформулировал совершенно "по другому", исказив смысл исходной задачи
Вот этот "тезис":

Цитата vitaly333 @ 03.10.09, 11:51

Если в xmm0 записать {A1,A1} а в xmm1 - {B1,B2}

вообще не имеет отношения к исходной задаче, т.к. ничего не дает кроме лишних затрат процессорного времени. SSE операции особенно эффективны, когда входные данные ипользуют "как есть" без всяких перестановок. Применительно к данному случаю - если даны строка A=(A1,A2,..) и строка B=(B1,B2,..), то для получения суммы произведений A[k]*B[k] нужно вычислять (С',C")= (A[k],A[k+1])*(B[k],B[k+1]); k++=2, где C' и C" - это частичные суммы по четным и нечетным элементам, которые складываются друг с другом только по окончании цикла. Для повышения эффективности можно в одном цикле умножать одну строку A на 2 или 4 строки B - именно этот вариант я и привел в #6. (В принципе можно и наоборот: умножать одну B на на 2 или 4 A, но при этом выходные C будут идти не по порядку - потенциально хуже кэширование).

Цитата vitaly333 @ 04.10.09, 09:30

Имеет.Если вы под "перемножением в лоб" подрузамеваете алгортим строка-строка, то он работает эффективно только до тех пор пока 3 строки умещаются к кэш. Блочный же алгортим эффективен на любой размерности матриц A,B,C.

Не забывай, что кэш это маленький сверхбыстрый L1 + достаточно большой, но не такой шустрый L2. Ты почему-то "бездумно" зациклился именно на L1 и удивляешься "почему то не работает префетч". А он вообще должен работать в данных условиях ? Префетчить блок C вообще не имеет никакого смысла, т.к. каждый элемент C вычисляется сравнительно долго и соотв-но линейки блока C могут безболезненно подгружаться даже из ОЗУ, не говоря уже о L2. Элементы блоков A и B используются чаще, но в любом случае для просчета всего блока C требуется по многу раз пройтись во всем элментам блоков A и B, поэтому если и будут тормоза на первом проходе, то они нивелируются при последующих проходах, когда данные уже сами по себе будут в кэше.
Вывод: про L1 в данном случае можно забыть и оптимизировать блоки только в расчете на L2, а это уже более чем на порядок бОльшие размеры. По крайней мере для твоего примера c N ~ 2000 в L2 запросто умещаются десятки строк. Тем не менее блочная обработка видимо по любому нужна, но скорее всего не "квадратно-гнездовая" с одинаковыми размерами по C и (A,B)

Цитата vitaly333 @ 04.10.09, 09:30

Как получу ответ на 1 вопрос

Я уже предложил вариант #6 (осталось только объединение половинок прикрутить) и обосновал его чуть выше. Ждемс реакции

Сообщение отредактировано: leo - 04.10.09, 10:53

[leo]

Можно и твой вариант использовать с попарным умножением двух строк A на две строки B, только в "нормальной" SSE-реализации он должен выглядеть примерно так:

  ;eax = *A1; edx = *A2; edi = *B1; esi = *B2;
  ;ecx = N*8; eax=eax+ecx; ... esi=esi+ecx
  ;neg ecx
  ;xmm0 = xmm1 = xmm2 = xmm3 = 0 - накапливаемые суммы для C1..C4
loop1:
  movapd xmm4,[eax+ecx]  ;A1
  movapd xmm5,xmm4
  movapd xmm6,[edi+ecx]  ;B1
  movapd xmm7,[esi+ecx]  ;B2
  mulpd  xmm4,xmm6       ;A1*B1
  mulpd  xmm5,xmm7       ;A1*B2
  addpd xmm0,xmm4
  addpd xmm1,xmm5
  movapd xmm4,[edx+ecx]  ;A2
  mulpd  xmm6,xmm4       ;A2*B1
  mulpd  xmm7,xmm4       ;A2*B2
  addpd xmm2,xmm6
  addpd xmm3,xmm7
  add ecx,16
  jnz loop1

Добавлено 05.10.09, 05:32
Ну и объединение результатов после завершения цикла:

  mov ecx,pC1       ;ecx = *C1
  movapd xmm4,xmm0
  unpklpd xmm0,xmm1 ;(C1',C2')
  unpkhpd xmm4,xmm1 ;(C1",C2")
  addpd xmm4,xmm0   ;(C1,C2)
  addpd [ecx],xmm5
  ;---------------
  mov ecx,pC3  ;ecx = *C3 для варианта #11 или add ecx,16 для варианта #6
  movapd xmm5,xmm2
  unpklpd xmm2,xmm3 ;(C3',C4')
  unpkhpd xmm5,xmm3 ;(C3",C4")
  addpd xmm5,xmm2   ;(C3,C4)
  addpd [ecx],xmm5

Сообщение отредактировано: leo - 05.10.09, 05:32

[vitaly333]

leo
Попробовал свой вариант и ваши оба варианта SSE2 векторизации с умножением 1 строки матрицы A на 4 строки матрицы B (ваш пост №6) и с умножением 2 строк A на две строки матрицы B(пост №11). Мой вариант "хлам". №11 несколько быстрее №6. Я так полагаю за счет меньшего количества загрузок из кэша. Впринципе уже не плохо.
Для ещё большего ускорения алгоритма можно развернуть вручную внутренний цикл на n-ое кол-во раз чтобы избежать тормозов конвеера. Как определить оптимальное количество разворотов и как это вообще сделать технически(для каждого кол-ва разворотов писать фун-цию на ASM-е и потом инлайнить её в вызывающую ф-ю)?

Цитата

Не забывай, что кэш это маленький сверхбыстрый L1 + достаточно большой, но не такой шустрый L2. Ты почему-то "бездумно" зациклился именно на L1 и удивляешься "почему то не работает префетч". А он вообще должен работать в данных условиях ? Префетчить блок C вообще не имеет никакого смысла, т.к. каждый элемент C вычисляется сравнительно долго и соотв-но линейки блока C могут безболезненно подгружаться даже из ОЗУ, не говоря уже о L2. Элементы блоков A и B используются чаще, но в любом случае для просчета всего блока C требуется по многу раз пройтись во всем элментам блоков A и B, поэтому если и будут тормоза на первом проходе, то они нивелируются при последующих проходах, когда данные уже сами по себе будут в кэше.

А в каких же тогда случаях полезен префетч? Как его можно применить к задаче умножения матриц?

Цитата

Вывод: про L1 в данном случае можно забыть и оптимизировать блоки только в расчете на L2, а это уже более чем на порядок бОльшие размеры. По крайней мере для твоего примера c N ~ 2000 в L2 запросто умещаются десятки строк. Тем не менее блочная обработка видимо по любому нужна, но скорее всего не "квадратно-гнездовая" с одинаковыми размерами по C и (A,B)

Т.е. размер блока нужно выбирать таким чтобы 3 блока умещались в L2? Знаете, я провел небольшой эксперимент. Для задачи N=4032 я взял сначала размер блока=48 (3 блока умещаются в мой L1), 96(где то между L1 и L2)и потом 144(3 блока помещаются в L2). Результаты честно говоря не сильно различаются.

размер матриц = 4032x4032
Компилятор Intel C++ v11.1 /O3:
Метод: blockMult(без SSE)
blockSize = 48 Время = 63.3
blockSize = 96 Время = 58.0
blockSize = 144 Время = 59.0

Метод: blockMult(с ручным SSE)
blockSize = 48 Время = 51.7
blockSize = 96 Время = 50.0
blockSize = 144 Время = 52.2

Компилятор GNU G++ v4.4 -O3:
Метод: blockMult(без SSE)
blockSize = 48 Время = 68.7
blockSize = 96 Время = 68.4
blockSize = 144 Время = 69

Так что как выбирать оптимальный размер блока для меня пока загадка.

Сообщение отредактировано: vitaly333 - 08.10.09, 20:09

[leo]

Цитата vitaly333 @ 08.10.09, 20:05

Для ещё большего ускорения алгоритма можно развернуть вручную внутренний цикл на n-ое кол-во раз чтобы избежать тормозов конвеера

Тут уже и так все нормально развернуто, дальше разворачивать не имеет смысла. Кстати учти, что 128-битные packed операции только в Core2 и K10 выполняются за 1 микрооперацию, а во всех предшествующих процах за 2 мопа (половинками по 64 бита). Поэтому тут и так получается около 25 мопов SSE2 на цикл - больше некуда. Параллельного вычисления 4-х сумм также достаточно, чтобы компенсировать ("скрыть") латентность addpd

Цитата vitaly333 @ 08.10.09, 20:05

А в каких же тогда случаях полезен префетч?

В специфических
Классический пример - копирование больших блоков памяти с использованием movntX. В этом случае блок-префетч улучшает работу шины памяти и ОЗУ, т.к. данные сначала только читаются из ОЗУ, а затем только пишутся. В противном случае идет чередование чтения\записи и возникают доп.тормоза

Цитата vitaly333 @ 08.10.09, 20:05

я провел небольшой эксперимент. Для задачи N=4032 я взял сначала размер блока=48 (3 блока умещаются в мой L1), 96(где то между L1 и L2)и потом 144(3 блока помещаются в L2). Результаты честно говоря не сильно различаются

Я тебе уже дважды намекал, что брать большие квадратные блоки A и B это далеко не оптимально. Ты не учитываешь хардварный префетч (ХВП) и увеличением числа просчитываемых строк нарушаешь его работу. Не знаю сколько одновременных потоков ХВП поддерживают атлоны (нет инфы), а вот у Intel все расписано - до 8 потоков в P4 и до 16 в Core. А ты со своим 48х48 блоком пытаешься одновременно работать с 48*2 = 96 строками, т.е. с 96 разными потоками данных одновременно - ес-но ХВП от этого дуреет и отдыхает на лавочке
Блоки A и B нужно брать не квадоратные, а "строчные", т.е. во всю длину строки и просчет вести по всей длине строки - при этом и ХВП будет нормально работать и своевремменно подгружать данные как из ОЗУ в L2, так и из L2 в L1. Прогнал один полный просчет 2х2 - получил 4 готовых значения С, переходишь к следующей итерации - оставляешь те же 2 строки A и 2 новые строки B (чтобы двигаться по строкам C для просчета полной линейки кэша). Число просчитываемых строк B выбираем из условия, чтобы в L2 умещалось одинаковое число строк A и B. Например, при размере N = 4032 получаем 8 строк для L2=512К или 16 для 1Mb (как раз на одну или 2 линейки С). Просчитали нужное число строк по B - берем следующие 2 строки A и все повторяем. В итоге когда просчитали все 8 или 16 строк по А получаем законченный квадратный блок C размером 8х8 или 16х16. Переходим к следующему 8-строчному блоку по B с тем же блоком А.
В итоге при таком построчном просчете данные подгружаются последовательно и не более чем 4-мя потоками ХВП одновременно, а бОльшую часть времени всего 2 потоками - всё дубово-линейно и все довольны

[leo]

Цитата leo @ 09.10.09, 14:27

В итоге когда просчитали все 8 или 16 строк по А получаем законченный квадратный блок C размером 8х8 или 16х16. Переходим к следующему 8-строчному блоку по B с тем же блоком А.

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

PS: Тут видимо нужно звать на помощь спецов-алгоритмистов (albom, отзовись, а-у-у )
У меня лично моСК отказывается варить, и подло намекает, что без непрерывной подгрузки одной из матриц из ОЗУ эта задачка не решается и соотв-но надежды на какую-то оптимальную блочную обработку - призрачны и тщетны

Сообщение отредактировано: leo - 09.10.09, 19:57

[leo]

Цитата leo @ 09.10.09, 19:56

без непрерывной подгрузки одной из матриц из ОЗУ эта задачка не решается

Тут я ес-но опять не прав. Если грузить блок B и затем умножать его на все строки A, то загрузка из памяти происходит только при
1) переходе к новому блоку B (грузится весь блок)
2) переходе к новой паре строк A (грузится пара строк A при умножении на первую пару строк B)
3) редко и по мелочи при сохранении линеек C (2 линейки на 2 строки A и 8 строк B)
Т.е. во время умножения пары строк A на вторую и последуюшие пары строк B работа идет только с кэшем (не считая сохранения линеек С). Поэтому чем больше блок B, тем лучше, т.е. его нужно брать почти на весь размер L2 (оставив место под 2 строки A и несколько линеек C). При этом лучше брать число строк, кратное 8 (в последнем - как получится), чтобы получать на выходе целое число линеек С.
Т.е. в итоге получается, что работа идет горизонтальными полосам по B и вертикальными по C, и для каждой такой полосы пробегается вся мартица A (т.е.блочности по A как таковой нет)

[vitaly333]

leo, я попробовал ваши варианты строчного блокирования. Оба они оказались хуже чем вариант с квадратными блоками одного размера.Вот результаты для той же задачи с N = 4032:
1) Размер блока = 8(Это когда 8 строк матрицы A умножаются на 8 строк матрицы B) Время = 89.4
2) Размер блока = 14 (Это когда полоска матрицы B умножаеся на всю матрицу A) Время = 81.3
Поясню как я вычислял размер блока. При N = 4032 каждая строка будет занимать 4032*8 = 32256 байт. Как вы писали, во 2-ом случае размер блока нужно брать почти на весь размер L2,оставив место под 2 строки A и несколько линеек C. При L2 = 512 000, память, необходимая для хранения блока B = L2 -2*32256 = 447488. Размер блока = 447478/32256 = 13.8
Значит дело здесь не в работе ХВП а в чем то другом. Я думаю стоит посмотреть в сторону оптимизированных библиотек BLAS, таких как Atlas,MKL,GOTOBlas и др. Какие они используют схемы блокирования и как вообще они там выбирают размеры блоков. На затравку вот нашел статью ведущего инженера Intel MKL Грега Генри. Там, в частности пишется что блокирование нужно осуществлять на всех уровнях иерархии памяти, на каждом переходе от менее скоростной памяти к более скоростной.

Цитата

А вообще, vitaly333, покажи тот вариант, который у тебя сейчас получился (если это будет чем-то вроде функции на C, которую можно сразу запустить и замерить скорость, — просто замечательно), потому как сообщение о том, что матрица 4032x4032 перемножается за 50.0 секунд (это секунды ведь?) может говорить либо о том, что у тебя плохой алгоритм/реализация, либо о том, что у тебя процессор старый

Незнаю, процессор вроде не такой уж и старый: AMD Athlon64x2 2.41 GHZ.
Вот вариант:

#include "stdafx.h"
#include "conio.h"
#include "math.h"
#include <iostream>
#include "time.h"
#include "malloc.h"
 
using namespace std;
 
/*
Ограничения:
- допустимы только квадратные матрицы
- один поток
Выбор размера блока:
- должен быть кратен N
- желательно чтобы был кратен 8 и был <=sqrt(L/24*bsize), где L - размер кэша L1 или L2 в байтах
*/
double** blockMultSSE2(double**A, double** B,double**C, int N,int bsize){
 
    int n,m,iter,i,j,k;
    int ibeg,iend,jbeg,jend,kbeg,kend;
    
     double *a1,*a2,*b1,*b2,*c1,*c2;
 
     transpose(B,N);
 
     int addrC = 0,addrBeg = 0;
     int gridSize = N/bsize;
 
 
    for (n=0; n<gridSize; n++){
     ibeg = n*bsize;
     iend = ibeg+bsize;
     for (m=0; m<gridSize; m++){
       jbeg = m*bsize;
       jend = jbeg+bsize;                  
        for (iter=0; iter<gridSize; iter++) {
        kbeg = iter*bsize;
        kend = kbeg+bsize;
        addrBeg = kbeg*8;
          for (i=ibeg; i<iend; i+=2){
            a1 = A[i];
            a2 = A[i+1];
            c1 = C[i];
            c2 = C[i+1];                            
             for (j=jbeg; j<jend; j+=2){
                b1 = B[j];
                b2 = B[j+1];
                addrC = j*8;
                 _asm{
                                   mov eax,a1
                           mov edx,a2
                                   mov edi,b1
                                   mov esi,b2
                                
                           mov ecx,bsize ; итератор цикла
                               mov ebx,addrBeg
                                
                                  ;xmm0 = xmm1 = xmm2 = xmm3 = 0 - накапливаемые суммы для C1..C4
                   xorps xmm0,xmm0
                   xorps xmm1,xmm1
                   xorps xmm2,xmm2
                   xorps xmm3,xmm3
                          
                                 loop1:
 
                                   movapd xmm4,[eax+ebx]  ;A1
                                   movapd xmm5,xmm4
                                   movapd xmm6,[edi+ebx]  ;B1
                                   movapd xmm7,[esi+ebx]  ;B2
                                   mulpd  xmm4,xmm6       ;A1*B1
                                   mulpd  xmm5,xmm7       ;A1*B2
                                   addpd xmm0,xmm4
                                   addpd xmm1,xmm5
                                   movapd xmm4,[edx+ebx]  ;A2
                                   mulpd  xmm6,xmm4       ;A2*B1
                                   mulpd  xmm7,xmm4       ;A2*B2
                                   addpd xmm2,xmm6
                                   addpd xmm3,xmm7
                   add ebx,16
                   sub ecx,2
 
                                 jnz loop1
                            
                   mov ebx,addrC
 
                                   mov ecx,c1       ;ecx = *C1                            
                                   movapd xmm4,xmm0                        
                   movapd xmm6,[ecx+ebx]
                                   unpcklpd xmm0,xmm1
                                   unpckhpd xmm4,xmm1
                                   addpd xmm4,xmm0  
                                   addpd xmm6,xmm4
                                   movapd [ecx+ebx],xmm6
  
                                   mov ecx,c2        ;ecx = *C2
                                   movapd xmm5,xmm2
                   movapd xmm6,[ecx+ebx]
                                   unpcklpd xmm2,xmm3
                                   unpckhpd xmm5,xmm3
                                   addpd xmm5,xmm2  
                                   addpd xmm6,xmm5
                                   movapd [ecx+ebx],xmm6
                              
                    
                }
 
               }
            }
         }
       }
  }
    return C;
  }
 
void fillMatrices(double** A,double** B,double** C, int N){
 
        int i,j,k;
        k = 10.0;
         for (i = 0;i < N; i++) {
            for (j = 0;j<N;j++){
                A[i][j] = rand()*k;
                B[i][j] = rand()*k;
                C[i][j] = 0.0;
            }
         }
}
 
void transpose(double** A,int N){
 
int i,j;
double tmp;
 
 for (i=0;i<N;i++){
   for (j=i;j<N;j++){
     tmp = A[i][j];
     A[i][j] = A[j][i];
     A[j][i] = tmp;
   }
 }
 
}
 
 
int main() {
    
  const int N = 4096;
  const int BSIZE = 64;
 
  std::cout <<"Allocate memory..."<<endl;
 
  double** A = new double*[N];
  double** B = new double*[N];
  double** C = new double*[N];
  for (int i=0;i<N;i++){
    A[i] = (TYPE *)_aligned_malloc(N*sizeof(double),16);
    B[i] = (TYPE *)_aligned_malloc(N*sizeof(double),16);
    C[i] = (TYPE *)_aligned_malloc(N*sizeof(double),16);
  }
  std::cout <<"Fill matrices..."<<endl;
    
  fillMatrices(A,B,C,N);
 
  std::cout <<"Matrix multiplication..."<<endl;
 
  float t1,t2;
  t1 = clock();
  blockMult(A,B,C,N,BSIZE);
  t2 = (float)(clock()-t1)/CLOCKS_PER_SEC;
 
  std::cout <<"Computation time = "<<t2<<" sec. "<<endl;
  if (t2==0.0) t2 = 0.000000000001;
  float perform = 2*powf(N/100,3)/t2;
  std::cout<<"Performance = "<<perform<<" MFlops "<<endl;
  //getch();
 return 0;
}

albom, попробовал я и вашу реализацию. Ваш метод требует очень много памяти. С установленным на моем компьютере 1 GB мне не удалось посчитать задачку N = 4032. Винда попала в глубочайший своп. Поэтому привожу результаты для N = 2048 Intel Compiler v 11.1.048 -O3:
dmm - 7.1 сек.
blockMultSSE2 - 6.6 сек.

Интересно бы было посмотреть какие у вас получатся результаты.

Добавлено 14.10.09, 16:35
leo, я попробовал ваши варианты строчного блокирования. Оба они оказались хуже чем вариант с квадратными блоками одного размера.Вот результаты для той же задачи с N = 4032:
1) Размер блока = 8(Это когда 8 строк матрицы A умножаются на 8 строк матрицы B) Время = 89.4
2) Размер блока = 14 (Это когда полоска матрицы B умножаеся на всю матрицу A) Время = 81.3
Поясню как я вычислял размер блока. При N = 4032 каждая строка будет занимать 4032*8 = 32256 байт. Как вы писали, во 2-ом случае размер блока нужно брать почти на весь размер L2,оставив место под 2 строки A и несколько линеек C. При L2 = 512 000, память, необходимая для хранения блока B = L2 -2*32256 = 447488. Размер блока = 447478/32256 = 13.8
Значит дело здесь не в работе ХВП а в чем то другом. Я думаю стоит посмотреть в сторону оптимизированных библиотек BLAS, таких как Atlas,MKL,GOTOBlas и др. Какие они используют схемы блокирования и как вообще они там выбирают размеры блоков. На затравку вот нашел статью ведущего инженера Intel MKL Грега Генри. Там, в частности пишется что блокирование нужно осуществлять на всех уровнях иерархии памяти, на каждом переходе от менее скоростной памяти к более скоростной.

[vitaly333]

Что-то какие-то проблемы с форумом...

Сообщение отредактировано: vitaly333 - 14.10.09, 18:00