Само умножение нехитрое, элементы строк умножаются на столбцы поэлементно и складываются. Как корректно умножать можно посмотреть здесь Языковая часть написана на 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.
Возможно данная заметка поможет сэкономить время. Удачи!
Автор: Еремин Дмитрий
