Быстрое умножение матриц
На днях в блоге Zealint был объявлен конкурс на быстрое умножение матриц. Мне показалось это интересным и я принял в нем участие. Приведу основные моменты из своего кода. Может быть кому-то это пригодится.

Условия: Используются квадратные матрицы. Элементы матриц - целые числа, по модулю не превышающие 600. Максимальный порядок матрицы n = 5000. Исходные матрицы читаются из файла. Результат также выводится в файл.

Для эффективного использования кэша процессора матрица разбивается на матрицы меньшего размера. В [1] приводится вариант разбиения с изменением индекса. У меня использован другой вариант. Каждая меньшая матрица физически расположена в непрерывном участке памяти и все три матрицы представлены как матрицы, состоящие из M2xM2 матриц меньшего размера M1xM1:

pascal
type
  tm1 = array [0..M1-1, 0..M1-1] of smallint;
  tm2 = array [0..M2-1, 0..M2-1] of tm1;

Это объявление типа для исходных матриц A и B. Для матрицы результата объявление такое-же, но тип элементов integer. Естественно, что теперь размер исходных матриц может принимать только ряд дискретных значений и с лишними частями нужно что-то делать. У меня они просто заполняются нулями и используются в вычислениях.

Что касается констант, то после нескольких экспериментов было выбрано M1 = 64 и M2 = 79 (Все три меньших матрицы помещаются в кэш первого уровня, а 64*79=5056, что покрывает максимальный размер матрицы). Такое разбиение не очень сильно усложняет ввод/вывод матриц.

Теперь о вычислениях. Математически формула умножения матриц ничем не отличается для элементов, которые сами являются матрицами. Отсюда код:

pascal
procedure MulM1xM1(var R:tmI1; var A, B:tm1);
var i, j, k:integer;
    s0 : integer;
begin
  for i := 0 to M1-1 do begin
     for j := 0 to M1-1 do begin
        s0 := 0;
        k := 0;
        while k < M1 do begin
           s0 := s0 + A[i,k+ 0]*B[j,k+ 0]+
                      // ....   ....
                      A[i,k+15]*B[j,k+15];
           inc(k, 16)
        end;
        R[i,j] := R[i,j] + s0;
     end;
  end;
end;

procedure MulMatrix;
var i, j, k:integer;
begin
  for i := 0 to M-1 do begin
     for j := 0 to M-1 do begin
        for k := 0 to M-1 do begin
           MulM1xM1(R[i, j], A[i, k], B[j, k]);
        end;
     end;
  end;
end;

Здесь $ M = \lceil N / M1 \rceil $ - количество маленьких матриц в строке большой матрицы для заданного размера матрицы $ N $

Обратите внимание, для матрицы B изменен порядок следования индексов и, соответственно, используется транспонированная матрица. (Транспонирование производится при чтении и не сказывается на общем времени вычислений). Такой порядок индексов позволяет читать элементы массива по строкам, что значительно эффективнее. Кроме того, самый внутренний цикл частично развернут. Это еще один из способов оптимизации за счет устранения точек ветвления.

comment
Замечу, что прямая реализация на ассемблере ничего не дает в плане ускорения работы. Delphi генерирует достаточно эффективный код. В частности, при обращении к элементам массива не используется умножение, а делается что-то что-то вроде:

asm
 movsx eax, [esi+ebx*2+0]
 movsx edx, [edi+ebx*2+0]
 imul  eax, edx
 ;...

здесь ebx - счетчик k внутреннего цикла, 2-ка - размер элемента (элементы матрицы хранятся как signed word).

Далее самый внутренний цикл (while k) заменен ассемблерной вставкой с использованием команд MMX. На самом деле была заменена вся процедура, но здесь я приведу только код внутреннего цикла:

pascal
   asm
      mov  esi, PA  // @A[i, 0]
      mov  edi, PB  // @B[j, 0]

      movdqa  xmm0, dqword ptr [esi+ 00]  
      movdqa  xmm1, dqword ptr [esi+ 16]  
      // и т.д. для xmm2-xmm6 .... для 32, 48, 64, 80, 96
      movdqa  xmm7, dqword ptr [esi+112]  

      PMADDWD xmm0, dqword ptr [edi+ 00]
      PMADDWD xmm1, dqword ptr [edi+ 16]
      // и т.д. для xmm2-xmm6 .... для 32, 48, 64, 80, 96
      PMADDWD xmm7, dqword ptr [edi+112]

      paddd   xmm0, xmm1
      paddd   xmm0, xmm2
      // .... для xmm3 .. xmm7 
      paddd   xmm0, xmm7


      pshufd  xmm2, xmm0,1110b 
      paddd   xmm0, xmm2       
      pshufd  xmm2, xmm0,0001b 
      paddd   xmm0, xmm2

      movd    s0, xmm0
   end;

Как видим этот блок вычисляет произведение одной строки матрицы 64x64 на строку другой матрицы 64x64.

comment
При реализации столкнулся с тем, что в Delphi нет директивы выравнивания по границе параграфа. Максимум Align8. Пришлось выравнивать вручную добавлением лишних переменных.

Последние команды - сложение по горизонтали можно заменить командами PHADDD:

pascal
    // Сложение по горизонтали
    pxor    xmm1, xmm1
    //PHADDD xmm0, xmm1
    db $66, $0F, $38, $02, $C1
    //PHADDD xmm0, xmm1
    db $66, $0F, $38, $02, $C1

comment
К сожалению, Delphi о них не знает, да и время выполнения почти не меняется. Кстати, эти команды отказался выполнять имеющийся у меня (не самый старый) процессор AMD.

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

pascal
 asm
    PADDD  xmm0, dqword ptr [ebx]
    movdqa dqword ptr [ebx], xmm0
 end

Не буду приводить полный код процедуры, т.к. он достаточно большой, а основные части его уже приводились.


На этом все. Спасибо Артему (ака Zealint) за проделанную работу и мои поздравления venco с победой в конкурсе. Было интересно. Комментарии приветствуются.

1Intel 64 and IA-32 Architectures Optimization Reference Manual. Nov. 2007

20 февраля 2010

Comments
22.12.2010 01:18 mypepeemi
 Материал на пять с плюсом. Но есть и минус! У меня скорость интернета 56кб/сек. Страница грузилась около 40 секунд.

 
 
10.02.2016 02:34 shopwme.ru
 Вопрос о построении наиболее быстрого и устойчивого практического алгоритма 
 умножения больших матриц также остается нерешённым.
 
Вы можете оставить комментарий или задать вопрос
Ваше имя:

Текст сообщения:



Copyright © 2009-2014 by