Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Алгоритмы > Алгоритмы умножения матриц

Автор: Anark1 12.9.2008, 19:11

Интересуют подробные описания алгоритмов перемножения матриц.
Все, кроме классического правила (строка-столбец), алгоритма Штрассена, алгоритма Копперсмита-Винограда, если таковые еще существуют.

Автор: ksili 13.9.2008, 05:04

Не так давно проф. Грановский изобрёл алгоритм, который работает ещё быстрее Штрассена и Винограда. Сам алгоритм не знаю, видел его упоминание на форуме Intel. Интеловцы вроде его уже реализовали в своей Intel MKL.
В общем сейчас разработка новых алгоритмов идёт в основном в направлении распараллеливания. Алгоритм Грановского вроде тоже хорош для распараллеливания.

Ещё посмотри книги Воеводина. Он хорошо разбирает матрицы с т.зр. алгоритмов. Например, Воеводин, Кузнецов "Матрицы и вычисления".

Автор: Anark1 13.9.2008, 11:09

Описания алгоритма Грановского я не нашел, как впрочем и книгу... Она 1984 года, так что достать её весьма сложно.

Автор: ksili 13.9.2008, 11:25

Ну значит надо в библиотеку идти. Впрочем, как и за Воеводиным.
А я что-то думал, что алгоритм Грановского поновее будет, 90-х или даже 2000-х годов...

Автор: Anark1 13.9.2008, 12:09

Ладно, лови плюс, но тема все равно остается актуальной. smile

Если кто подскажет серьезные ресурсы по данной тематике, буду очень благодарен.

Автор: ksili 15.9.2008, 09:55

Вот http://num-anal.srcc.msu.ru/lib_na/libnal.htm можно посмотреть программы на Си и Фортране для различных алгоритмов умножения. В основных разделах выбирай "Линейная алгебра". Там будет пункт "Умножение матриц".

Автор: Anark1 15.9.2008, 20:41

ksili, мне это все не нужно. Программу я и так могу написать. Пока прицеливаюсь и расширяю кругозор smile

интересуют чисто алгоритмы, а там, насколько я понял, для вектор-матриц и прочих частных случаев используется частный случай алгоритма строка-столбец.

Автор: Vitaly333 15.10.2008, 22:57

Очень интересное наблюдение:

Обычный код для умножения двух квадратных матриц:

Код


for(i=0;i<N;i++)   {
     for(j=0;j<N;j++) {
         C(i,j)=0;
         for(k=0;k<N;k++) {
              C(i,j)+=A(i,k)*B(k,j);
          }
      }
  }

А теперь перед тем как умножить транспонируем матрицу B и изменим индексы:

Код


for(i=0;i<N;i++)   {
     for(j=0;j<N;j++) {
         C(i,j)=0;
         for(k=0;k<N;k++) {
              C(i,j)+=A(i,k)*B(j,k); 
          }
      }
  }

Это умножение будет в несколько раз (в 2-3 раза) быстрее чем стандартное!

Автор: ksili 16.10.2008, 04:16

1) Как замеряли время? (опишите эксперимент)
2) И почему там i большое и i маленькое?

Автор: Vitaly333 16.10.2008, 10:44

Цитата

1) Как замеряли время? (опишите эксперимент)

Замерял стандартными средствами языка!
Здесь всё дело в оптипальном расположении данных в кэше!

Цитата

2) И почему там i большое и i маленькое?

Очепятка вышла. Уже исправил!

Автор: ksili 16.10.2008, 14:23

Цитата(Vitaly333 @ 16.10.2008, 14:44

)

Здесь всё дело в оптипальном расположении данных в кэше!

Ну это уже будет зависеть от размера матриц.

К тому же мне кажется, что вы умножаете одним способом, а потом, в этой же программе, умножаете эти же матрицы. Т.е. ко второму умножению, обе матрицы оказываются в кэше! Понятно дело будет быстрее. Даже если оба умножения в разных прогах. То всё равно после транспонирования матрица В будет в кэше. Попробуйте замерить вместе с транспонированием.

Если ваша матрица B намного меньше кэша по размеру, то можно её без проблем поместить в кэш и без транспонирования. Хотя конечно последовательное чтение ячеек из области памяти, где находится матрица В, тоже играет свою роль.

Автор: DeNySkA 4.2.2009, 00:16

может кто знает или знает где искать как реализовать умножение 2 матриц на 3D структуре... Зарание спасибо!

Автор: bilbobagginz 8.2.2009, 00:54

Цитата(Vitaly333 @ 16.10.2008, 10:44

)

Замерял стандартными средствами языка!

а ответ будет равен в обеих "умножениях" ?
а то как-то A(i,k)*B(j,k) похоже на скалярное умножение двух строк, а не строка*столбик, как по определению...
поэкспериментируй с матрицами:
А=[1,2;3,4]; и B=[5,6;7,8]

Автор: Vitaly333 20.3.2009, 12:11

Цитата

Чем больше размер матриц тем больше будет видна разница в скорости. Оба умножения были запущены отдельно друг от друга и замеры времени я брал вместе с трансонированием. Дело в том, что матрица в современных системах программирования представляется как массив массивов. (Массив - столбец , каждый элемент которого, есть массив - строка). А массив данных - это какой то участок близко расположенных значений в памяти. Так вот процессору в матрице всегда "легче" пройтись по строчке, чем по столбцу, потому что, чтобы проходить по столбцу процессор каждый раз должен переключится на другой массив (другой участок в памяти) (как говорят обновить кэш) - а это дополнительное время. Поэтому в эффективных алгоритмах на матрицах нужно минимизировать такие переключения процессора.

Цитата

Конечно будет равен. Перед умножением матрица B - транспонируется.

Автор: W4FhLF 20.3.2009, 12:22

Да, я тоже такой приём использовал и у меня прирост на матрицах с количеством элементов большим 1 млн. составил ~1000%. На маленьких матрицах (< 256х256) на современных процессорах смысла в транспонировании нет, они целиком умещаются в кеш. Ещё при использовании такого приёма умножение можно векторизировать (т.е. заюзать SIMD и SSE), выполняя 2-4 умножения за такт.

А матрицу лучше хранить не строками, а линейно и адресовать как x[i * col + j]

Автор: Vitaly333 22.3.2009, 20:53

Цитата

Ещё при использовании такого приёма умножение можно векторизировать (т.е. заюзать SIMD и SSE), выполняя 2-4 умножения за такт.
А матрицу лучше хранить не строками, а линейно и адресовать как x[i * col + j]

Попробывал через векторизацию:

Код



int i,j,k,row,col,len;
        
    for (i = 0; i < n; i++) {
        row = i*n;
        len = row+n;
        for (j = row; j < len; j++) {    
            c[j] = 0.0;
            col = j-row;
            for (k = row; k < len; k++) {                
                c[j] += a[k] * b[col];
                col+=n;
            }
        }
    }

Но работает медленнее, чем 2-ой алгоритм с предварительным трансонированием.