Полный привод в матрицах 4×4

в 8:59, , рубрики: ассемблер, дельфи, математика

Само умножение нехитрое, элементы строк умножаются на столбцы поэлементно и складываются. Как корректно умножать можно посмотреть здесь Языковая часть написана на Delphi, а для оптимизации код выполнен с применением встроенного 64-х битного ассемблера. Рассматриваются 4 практических способа умножения.

Описание типа TMatrix4:

type TVec4 = array[0..3] of Single;
type TMatrix4 = array[0..3] of TVec4;

Получается матрица расположенная строками (ROW-major) и линейная в памяти — от 0-го элемента до 15-го. Из-за такой типизации выборка из массива на языке выглядит как Матрица[Ряд, столбец], т.е. первым стоит Y (ряд), а вторым X (столбец). Тип TVec4 применяется автором для извлечения отдельных строк (например, позиция = Матрица[3]) и не требует обязательного применения.

Скорость работы функций и процедур указаны в сводной таблице в конце.

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

Вариант A

function MulMatrixA(A, B: TMatrix4): TMatrix4;
var
  i, j: Integer;

begin
  for j := 0 to 3 do //Столбцы
  for i := 0 to 3 do //Строки
  begin
    Result[j, i] := A[j, 0] * B[0, i] +
                    A[j, 1] * B[1, i] +
                    A[j, 2] * B[2, i] +
                    A[j, 3] * B[3, i];
  end;
end;


Развернутый вариант выглядит сложно и запутанно, но выполняется быстрее на 40%.

Вариант B

function MulMatrixB(A, B: TMatrix4): TMatrix4;
var
  R: TMatrix4 absolute Result;

begin
  //Vec4 X
  R[0,0] := A[0,0] * B[0,0] + A[0,1] * B[1,0] + A[0,2] * B[2,0] + A[0,3] * B[3,0];
  R[0,1] := A[0,0] * B[0,1] + A[0,1] * B[1,1] + A[0,2] * B[2,1] + A[0,3] * B[3,1];
  R[0,2] := A[0,0] * B[0,2] + A[0,1] * B[1,2] + A[0,2] * B[2,2] + A[0,3] * B[3,2];
  R[0,3] := A[0,0] * B[0,3] + A[0,1] * B[1,3] + A[0,2] * B[2,3] + A[0,3] * B[3,3];

  //Vec4 Y
  R[1,0] := A[1,0] * B[0,0] + A[1,1] * B[1,0] + A[1,2] * B[2,0] + A[1,3] * B[3,0];
  R[1,1] := A[1,0] * B[0,1] + A[1,1] * B[1,1] + A[1,2] * B[2,1] + A[1,3] * B[3,1];
  R[1,2] := A[1,0] * B[0,2] + A[1,1] * B[1,2] + A[1,2] * B[2,2] + A[1,3] * B[3,2];
  R[1,3] := A[1,0] * B[0,3] + A[1,1] * B[1,3] + A[1,2] * B[2,3] + A[1,3] * B[3,3];

  //Vec4 Z
  R[2,0] := A[2,0] * B[0,0] + A[2,1] * B[1,0] + A[2,2] * B[2,0] + A[2,3] * B[3,0];
  R[2,1] := A[2,0] * B[0,1] + A[2,1] * B[1,1] + A[2,2] * B[2,1] + A[2,3] * B[3,1];
  R[2,2] := A[2,0] * B[0,2] + A[2,1] * B[1,2] + A[2,2] * B[2,2] + A[2,3] * B[3,2];
  R[2,3] := A[2,0] * B[0,3] + A[2,1] * B[1,3] + A[2,2] * B[2,3] + A[2,3] * B[3,3];

  //Vec4 W
  R[3,0] := A[3,0] * B[0,0] + A[3,1] * B[1,0] + A[3,2] * B[2,0] + A[3,3] * B[3,0];
  R[3,1] := A[3,0] * B[0,1] + A[3,1] * B[1,1] + A[3,2] * B[2,1] + A[3,3] * B[3,1];
  R[3,2] := A[3,0] * B[0,2] + A[3,1] * B[1,2] + A[3,2] * B[2,2] + A[3,3] * B[3,2];
  R[3,3] := A[3,0] * B[0,3] + A[3,1] * B[1,3] + A[3,2] * B[2,3] + A[3,3] * B[3,3];
end;

Оптимизированный вариант на ассемблере. Очень быстрый. Сделан для обычных переменных без выравнивания в памяти.

Вариант C

function MulMatrixC(A, B: TMatrix4): TMatrix4;
asm
  .NOFRAME

  //--------------------------------------------
  // Загружаем матрицу A строками
  //--------------------------------------------

  movups xmm0, dqword ptr [A]       // A[00..03]
  movups xmm1, dqword ptr [A + 16]  // A[04..07]
  movups xmm2, dqword ptr [A + 32]  // A[08..11]
  movups xmm3, dqword ptr [A + 48]  // A[12..15]

  //--------------------------------------------
  // Загружаем матрицу B строками
  //--------------------------------------------

  movups xmm8, dqword ptr [B]       // B[00..03]
  movups xmm9, dqword ptr [B + 16]  // B[04..07]
  movups xmm10, dqword ptr [B + 32] // B[08..11]
  movups xmm11, dqword ptr [B + 48] // B[12..15]

  //--------------------------------------------
  // Транспонируем матрицу B
  //--------------------------------------------

  // Прямой вид          | Вид в процессоре
  // [00] [01] [02] [03] | [03] [02] [01] [00]
  // [04] [05] [06] [07] | [07] [06] [05] [04]
  // [08] [09] [10] [11] | [11] [10] [09] [08]
  // [12] [13] [14] [15] | [15] [14] [13] [12]

  // Дублируем значения [00, 05, 10, 15]
  // Размер команды MOVAPS - 4 байта
  // что позволяет максимально быстро
  // копировать между регистрами
  // Для сравнения: MOVSD - 5 байтов

  movaps xmm4, xmm8  // [00]
  movaps xmm5, xmm9  // [05]
  movaps xmm6, xmm10 // [10]
  movaps xmm7, xmm11 // [15]

  // Первый столбец
  insertps xmm4, xmm9, 00010000b   // [04 -> 01]
  insertps xmm4, xmm10, 00100000b  // [08 -> 02]
  insertps xmm4, xmm11, 00110000b  // [12 -> 03]

  // Второй столбец
  insertps xmm5, xmm8, 01000000b   // [01 -> 04]
  insertps xmm5, xmm10, 01100000b  // [09 -> 06]
  insertps xmm5, xmm11, 01110000b  // [13 -> 07]

  // Третий столбец
  insertps xmm6, xmm8, 10000000b   // [02 -> 08]
  insertps xmm6, xmm9, 10010000b   // [06 -> 09]
  insertps xmm6, xmm11, 10110000b  // [14 -> 11]

  // Четвертый столбец
  insertps xmm7, xmm8, 11000000b   // [03 -> 12]
  insertps xmm7, xmm9, 11010000b   // [07 -> 13]
  insertps xmm7, xmm10, 11100000b  // [11 -> 14]

  //--------------------------------------------
  // Умножение матриц AxB
  //--------------------------------------------

  //------------------------
  // Матрица[0]
  //------------------------

  // Дублируем первую строку
  movaps xmm9, xmm0
  movaps xmm10, xmm0
  movaps xmm11, xmm0

  // Умножаем первую строку A на столбцы B[0..3]
  mulps xmm0, xmm4
  mulps xmm9, xmm5
  mulps xmm10, xmm6
  mulps xmm11, xmm7

  // Складываем все значения по горизонтали
  haddps xmm0, xmm0
  haddps xmm0, xmm0
  haddps xmm9, xmm9
  haddps xmm9, xmm9
  haddps xmm10, xmm10
  haddps xmm10, xmm10
  haddps xmm11, xmm11
  haddps xmm11, xmm11

  // Извлекаем результаты умножения
  movlhps xmm0, xmm10
  insertps xmm0, xmm9, 01010000b
  insertps xmm0, xmm11, 11110000b

  //------------------------
  // Матрица[1]
  //------------------------

  // Дублируем вторую строку
  movaps xmm9, xmm1
  movaps xmm10, xmm1
  movaps xmm11, xmm1

  // Умножаем вторую строку A на столбцы B[0..3]
  mulps xmm1, xmm4
  mulps xmm9, xmm5
  mulps xmm10, xmm6
  mulps xmm11, xmm7

  // Складываем все значения по горизонтали
  haddps xmm1, xmm1
  haddps xmm1, xmm1
  haddps xmm9, xmm9
  haddps xmm9, xmm9
  haddps xmm10, xmm10
  haddps xmm10, xmm10
  haddps xmm11, xmm11
  haddps xmm11, xmm11

  // Извлекаем результаты умножения
  movlhps xmm1, xmm10
  insertps xmm1, xmm9, 01010000b
  insertps xmm1, xmm11, 11110000b

  //------------------------
  // Матрица[2]
  //------------------------

  // Дублируем третью строку
  movaps xmm9, xmm2
  movaps xmm10, xmm2
  movaps xmm11, xmm2

  // Умножаем третью строку A на столбцы B[0..3]
  mulps xmm2, xmm4
  mulps xmm9, xmm5
  mulps xmm10, xmm6
  mulps xmm11, xmm7

  // Складываем все значения по горизонтали
  haddps xmm2, xmm2
  haddps xmm2, xmm2
  haddps xmm9, xmm9
  haddps xmm9, xmm9
  haddps xmm10, xmm10
  haddps xmm10, xmm10
  haddps xmm11, xmm11
  haddps xmm11, xmm11

  // Извлекаем результаты умножения
  movlhps xmm2, xmm10
  insertps xmm2, xmm9, 01010000b
  insertps xmm2, xmm11, 11110000b

  //------------------------
  // Матрица[3]
  //------------------------

  // Дублируем четвертую строку
  movaps xmm9, xmm3
  movaps xmm10, xmm3
  movaps xmm11, xmm3

  // Умножаем четвертую строку A на столбцы B[0..3]
  mulps xmm3, xmm4
  mulps xmm9, xmm5
  mulps xmm10, xmm6
  mulps xmm11, xmm7

  // Складываем все значения по горизонтали
  haddps xmm3, xmm3
  haddps xmm3, xmm3
  haddps xmm9, xmm9
  haddps xmm9, xmm9
  haddps xmm10, xmm10
  haddps xmm10, xmm10
  haddps xmm11, xmm11
  haddps xmm11, xmm11

  // Извлекаем результаты умножения
  movlhps xmm3, xmm10
  insertps xmm3, xmm9, 01010000b
  insertps xmm3, xmm11, 11110000b

  //------------------------
  // Результат
  //------------------------

  movups dqword ptr [Result], xmm0
  movups dqword ptr [Result + 16], xmm1
  movups dqword ptr [Result + 32], xmm2
  movups dqword ptr [Result + 48], xmm3
end;

Оптимизированный вариант на ассемблере. Самый быстрый. Сделан для переменных с выравниванием в памяти на 16 байт. И заметьте, это не функция, а процедура. Она возвращает результат в параметр RM, т.к. RM — это финальная переменная, а не внутренняя переменная Result, которая находится в стеке и из-за которой результат пишется дважды!

Доказательство про «пишется дважды»

Main.pas.1166: M^ := MulMatrixC(MC^, MD^);
000000000082DC36 488D8D60010000   lea rcx,[rsp+$00000160]
000000000082DC3D 488B9548020000   mov rdx,[rsp+$00000248]
000000000082DC44 4C8B8540020000   mov r8,[rsp+$00000240]
000000000082DC4B E820D7F1FF       call MulMatrixC //<-- Один раз пишется в функции MulMatrix C в [Result] т.е. стек
000000000082DC50 488BBD50020000   mov rdi,[rsp+$00000250]
000000000082DC57 488DB560010000   lea rsi,[rsp+$00000160]
000000000082DC5E C7C108000000     mov ecx,$00000008
000000000082DC64 F348A5           rep movsq //<-- Второй раз пишется после выхода из функции в переменную!

В MulMatrixD грузится указатель сразу на переменную, поэтому запись одна,
да и результат в таблице показателен.

Вариант D

procedure MulMatrixD(A, B, RM: TMatrix4);
asm
  .NOFRAME

  //--------------------------------------------
  // Загружаем матрицу A строками
  //--------------------------------------------

  movaps xmm0, dqword ptr [A]       // A[00..03]
  movaps xmm1, dqword ptr [A + 16]  // A[04..07]
  movaps xmm2, dqword ptr [A + 32]  // A[08..11]
  movaps xmm3, dqword ptr [A + 48]  // A[12..15]

  //--------------------------------------------
  // Загружаем матрицу B строками
  //--------------------------------------------

  movaps xmm8, dqword ptr [B]       // B[00..03]
  movaps xmm9, dqword ptr [B + 16]  // B[04..07]
  movaps xmm10, dqword ptr [B + 32] // B[08..11]
  movaps xmm11, dqword ptr [B + 48] // B[12..15]

  //--------------------------------------------
  // Транспонируем матрицу B
  //--------------------------------------------

  // Прямой вид          | Вид в процессоре
  // [00] [01] [02] [03] | [03] [02] [01] [00]
  // [04] [05] [06] [07] | [07] [06] [05] [04]
  // [08] [09] [10] [11] | [11] [10] [09] [08]
  // [12] [13] [14] [15] | [15] [14] [13] [12]

  // Дублируем значения [00, 05, 10, 15]
  // Размер команды MOVAPS - 4 байта
  // что позволяет максимально быстро
  // копировать между регистрами
  // Для сравнения: MOVSD - 5 байтов

  movaps xmm4, xmm8  // [00]
  movaps xmm5, xmm9  // [05]
  movaps xmm6, xmm10 // [10]
  movaps xmm7, xmm11 // [15]

  // Первый столбец
  insertps xmm4, xmm9, 00010000b   // [04 -> 01]
  insertps xmm4, xmm10, 00100000b  // [08 -> 02]
  insertps xmm4, xmm11, 00110000b  // [12 -> 03]

  // Второй столбец
  insertps xmm5, xmm8, 01000000b   // [01 -> 04]
  insertps xmm5, xmm10, 01100000b  // [09 -> 06]
  insertps xmm5, xmm11, 01110000b  // [13 -> 07]

  // Третий столбец
  insertps xmm6, xmm8, 10000000b   // [02 -> 08]
  insertps xmm6, xmm9, 10010000b   // [06 -> 09]
  insertps xmm6, xmm11, 10110000b  // [14 -> 11]

  // Четвертый столбец
  insertps xmm7, xmm8, 11000000b   // [03 -> 12]
  insertps xmm7, xmm9, 11010000b   // [07 -> 13]
  insertps xmm7, xmm10, 11100000b  // [11 -> 14]

  //--------------------------------------------
  // Умножение матриц AxB
  //--------------------------------------------

  //------------------------
  // Матрица[0]
  //------------------------

  // Дублируем первую строку
  movaps xmm9, xmm0
  movaps xmm10, xmm0
  movaps xmm11, xmm0

  // Умножаем первую строку A на столбцы B[0..3]
  mulps xmm0, xmm4
  mulps xmm9, xmm5
  mulps xmm10, xmm6
  mulps xmm11, xmm7

  // Складываем все значения по горизонтали
  haddps xmm0, xmm0
  haddps xmm0, xmm0
  haddps xmm9, xmm9
  haddps xmm9, xmm9
  haddps xmm10, xmm10
  haddps xmm10, xmm10
  haddps xmm11, xmm11
  haddps xmm11, xmm11

  // Извлекаем результаты умножения
  movlhps xmm0, xmm10
  insertps xmm0, xmm9, 01010000b
  insertps xmm0, xmm11, 11110000b

  //------------------------
  // Матрица[1]
  //------------------------

  // Дублируем вторую строку
  movaps xmm9, xmm1
  movaps xmm10, xmm1
  movaps xmm11, xmm1

  // Умножаем вторую строку A на столбцы B[0..3]
  mulps xmm1, xmm4
  mulps xmm9, xmm5
  mulps xmm10, xmm6
  mulps xmm11, xmm7

  // Складываем все значения по горизонтали
  haddps xmm1, xmm1
  haddps xmm1, xmm1
  haddps xmm9, xmm9
  haddps xmm9, xmm9
  haddps xmm10, xmm10
  haddps xmm10, xmm10
  haddps xmm11, xmm11
  haddps xmm11, xmm11

  // Извлекаем результаты умножения
  movlhps xmm1, xmm10
  insertps xmm1, xmm9, 01010000b
  insertps xmm1, xmm11, 11110000b

  //------------------------
  // Матрица[2]
  //------------------------

  // Дублируем третью строку
  movaps xmm9, xmm2
  movaps xmm10, xmm2
  movaps xmm11, xmm2

  // Умножаем третью строку A на столбцы B[0..3]
  mulps xmm2, xmm4
  mulps xmm9, xmm5
  mulps xmm10, xmm6
  mulps xmm11, xmm7

  // Складываем все значения по горизонтали
  haddps xmm2, xmm2
  haddps xmm2, xmm2
  haddps xmm9, xmm9
  haddps xmm9, xmm9
  haddps xmm10, xmm10
  haddps xmm10, xmm10
  haddps xmm11, xmm11
  haddps xmm11, xmm11

  // Извлекаем результаты умножения
  movlhps xmm2, xmm10
  insertps xmm2, xmm9, 01010000b
  insertps xmm2, xmm11, 11110000b

  //------------------------
  // Матрица[3]
  //------------------------

  // Дублируем четвертую строку
  movaps xmm9, xmm3
  movaps xmm10, xmm3
  movaps xmm11, xmm3

  // Умножаем четвертую строку A на столбцы B[0..3]
  mulps xmm3, xmm4
  mulps xmm9, xmm5
  mulps xmm10, xmm6
  mulps xmm11, xmm7

  // Складываем все значения по горизонтали
  haddps xmm3, xmm3
  haddps xmm3, xmm3
  haddps xmm9, xmm9
  haddps xmm9, xmm9
  haddps xmm10, xmm10
  haddps xmm10, xmm10
  haddps xmm11, xmm11
  haddps xmm11, xmm11

  // Извлекаем результаты умножения
  movlhps xmm3, xmm10
  insertps xmm3, xmm9, 01010000b
  insertps xmm3, xmm11, 11110000b

  //------------------------
  // Результат
  //------------------------

  movaps dqword ptr [RM], xmm0
  movaps dqword ptr [RM + 16], xmm1
  movaps dqword ptr [RM + 32], xmm2
  movaps dqword ptr [RM + 48], xmm3
end;

Итоговая таблица по вариантам:
(в скобках указан процент скорости от варианта D).

A – 7 376 698 умножений в секунду (16,4%)
B – 12 578 616 умножений в секунду (28,1%)
C – 37 735 849 умножений в секунду (84,3%)
D – 44 754 744 умножений в секунду (100%)

Получается вариант D быстрее A в 6 с лишним раз!

В некоторых случаях транспонирование матрицы B не требуется. Это значит, что блок транспонирования в версии D можно убрать, и тогда получается 5-й вариант:

E – 52 543 086 умножений в секунду (117,4%), что более чем в 7 раз выше скорости A!

Тестирование проводилось на процессоре Intel Core i7 6950X в одном потоке.
Компилятор Delphi XE6 Professional.

Возможно данная заметка поможет сэкономить время. Удачи!

Автор: Еремин Дмитрий

Источник


* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js