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

[vitaly333]

Цитата

Ну не знаю, от тебя ответы на почту дублируются и на форуме они некоторое врямя были видны (я же тебя цитировал). Может быть ты случайно у своих сообщений кнопку "Удалить" нажал раза 3-4 подряд?

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

Вот обновленные результаты.После обновления ваш метод стал значительно быстрее.

N	dmm2	my
1024X1024	0.755	0.797
2048X2048	5.532	6.468
4096X4096	39.161	64.210

Код своего блочного метода я прикрепил к письму. Так и не получается его сюда вставить.

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 - 14.10.09, 18:52

Прикреплённый файлBlockMult.rar (1.6 Кбайт, скачиваний: 112)

[vitaly333]

Цитата

но есть ведь более перспективный путь

Какой? Использовать алгоритмы типа Штрасcена?

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

[vitaly333]

Цитата

Да, использовать алгоритмы со сложностью ō(N3).

Можно попробывать, но это совсем другой путь. Хотелось бы сначала разобратся с блочными алгоритмами, ориентированными на кэш.

albomНе могли бы вы протестировать два варианта строчного блокирования, которые предлагал товарищ leo, а то после странных результов с Intel MKL я своей машине уже не доверяю. Код обоих методов в атаче. Для замеров времени можно взять ф-ию clock() для прошлого метода.

Сообщение отредактировано: vitaly333 - 15.10.09, 21:10

Прикреплённый файлrowBlockMults.rar (0.89 Кбайт, скачиваний: 109)

[leo]

Цитата vitaly333 @ 15.10.09, 21:09

Хотелось бы сначала разобратся с блочными алгоритмами, ориентированными на кэш.

Мне тоже. Причем разбираться желательно не путем "тупого" тестирования отдельных вариантов на конкретных тачках, а хоть с каким-то "теоретическим" обоснованием выбора размеров блока и способа прохода этих блоков (т.е. как раз организации "внешних циклов"). Что толку "тупо" сравнивать, к примеру, предложенный мною (на шару ) "построчный" вариант на атлонах с размером L2 512К и на Core2 c размером L2 в 2Мб или 4Мб - ежу понятно, что для такого способа обхода чем больше L2, тем лучше (на каждый блок требуется подгрузка из памяти двух строк матрицы A, поэтому чем больше строк в блоке, тем меньше относительный вклад подгрузки строк A в общее время расчета). К тому же на атлонах при той же частоте ОЗУ скорость чтения из памяти существенно ниже, чем в Core2
PS: Да и к тому же лобовой\тормозной вариант транспонирования матрицы может давать заметный вклад, поэтому нужно либо его привести к единому виду с оптимизированным вариантом albom'а, либо не учитывать транспонирование при тестировании

Добавлено 16.10.09, 05:30
albom
Можешь в двух словах пояснить суть твоего алгоритма dmm2 или придется разбираться по коду ?

[vitaly333]

Цитата

Можешь выложить сразу готовый к компиляции и запуску код/проект?

Могу.

Сообщение отредактировано: vitaly333 - 16.10.09, 08:31

Прикреплённый файлBlockMatrixMultiply.rar (4.54 Кбайт, скачиваний: 102)

[vitaly333]

Продолжение темы тут