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

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

Комментарии (39)


  1. picul
    13.10.2019 13:00

    Оптимизированный вариант — с insert'ами и hadd'ами? А с классическим вариантом сравнивали?


    1. eDmk Автор
      13.10.2019 14:27
      -1

      У мен нет C++. Не могу сравнить.


      1. picul
        13.10.2019 14:33

        Ну Вы можете и этот вариант реализовать в Delphi на ассемблере. Думаю, у Вас на это не уйдет много времени, а вот выигрыш в производительности скорее всего будет.


        1. eDmk Автор
          13.10.2019 14:36

          Это тоже самое что и у меня, только транспонирование идет за счет вычисления адресов.
          Т.е. загружается как транспонированная. У меня был такой вариант. Он работает немного медленнее. Перешел на транспонирование в регистрах.


          1. picul
            13.10.2019 14:55

            У Вас на каждые 4 умножения приходится по 8 hadd'ов (еще и неинтерливнутые для компенсации зависимостей, ай-яй-яй), весьма недешевых операций, между прочим. В варианте, который привел я, на каждые 4 умножения приходится по 3 обычных суммирования. В добавок еще и никаких транспонирований матриц не требуется (которые, кстати, спокойно выполняются за 8 shuffle'ов, не понимаю, зачем у Вас insert'ы).


            1. eDmk Автор
              13.10.2019 15:02

              В добавок еще и никаких транспонирований матриц не требуется (которые, кстати, спокойно выполняются за 8 shuffle'ов, не понимаю, зачем у Вас insert'ы).

              Это что то новое в умножении матриц :) Как можно зашуфлить вертикально?
              Транспонирование требуется всегда, за исключением случаев, когда матрица создается транспонированной.


              1. picul
                13.10.2019 15:21

                В C++ есть макрос "_MM_TRANSPOSE4_PS":

                Код
                #define _MM_TRANSPOSE4_PS(row0, row1, row2, row3) {                             __m128 _Tmp3, _Tmp2, _Tmp1, _Tmp0;                                                                                                          _Tmp0   = _mm_shuffle_ps((row0), (row1), 0x44);                      _Tmp2   = _mm_shuffle_ps((row0), (row1), 0xEE);                      _Tmp1   = _mm_shuffle_ps((row2), (row3), 0x44);                      _Tmp3   = _mm_shuffle_ps((row2), (row3), 0xEE);                                                                                          (row0) = _mm_shuffle_ps(_Tmp0, _Tmp1, 0x88);                          (row1) = _mm_shuffle_ps(_Tmp0, _Tmp1, 0xDD);                          (row2) = _mm_shuffle_ps(_Tmp2, _Tmp3, 0x88);                          (row3) = _mm_shuffle_ps(_Tmp2, _Tmp3, 0xDD);                      }


                1. eDmk Автор
                  13.10.2019 16:13

                  Я прекрасно знаю ассемблер и даже на опкодах могу писать.
                  При чем тут шуффлинг? Ваше умножение вообще не работает как умножение матриц. Вы нарушили весь порядок умножения и сложения.


                  1. picul
                    13.10.2019 16:55
                    +1

                    Я прекрасно знаю ассемблер и даже на опкодах могу писать.
                    Да я ж ничего не имею против)
                    При чем тут шуффлинг?
                    Я просто показал, как транспонируют матрицы с помощью shuffle'ов. Это лучше, чем транспонирование insert'ами.
                    Ваше умножение вообще не работает как умножение матриц.
                    Это не мое умножение, это общепринятый способ умножения, который используется в виду того, что он один из самых производительных из возможных в данных условиях. Тот код, который я привел — это не мой код, это одна из первых ссылок от гугления «sse 4x4 matrix multiplication». И не работает не он, а моя реализация на Delphi-ассемблере, которую я даже запустить не могу — в общем, не удивительно, что не работает.
                    Вы нарушили весь порядок умножения и сложения.
                    Какой еще порядок? Все так же, как у Вас, просто код быстрее.


                    1. eDmk Автор
                      13.10.2019 17:13
                      -1

                      Я просто показал, как транспонируют матрицы с помощью shuffle'ов. Это лучше, чем транспонирование insert'ами.

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


                      1. picul
                        13.10.2019 17:22
                        +1

                        Это транспонирование входит в стандартные заголовки всех распространенных компиляторов C++, и до Вас на него никто не жаловался. Перед тем, как называть его нерабочим, можно было бы скопировать его в Delphi и потестить, ну или разобраться, как оно работает.


                        1. eDmk Автор
                          13.10.2019 17:31
                          -1

                          Я вам свой вариант не навязывал. Я дал код для тех, кто пишет в Delphi, а C++ тут даже не упоминается. Я же не лезу с кодом на паскале в форум по Java? Вот и вы пришли совершенно не разбираясь в умножении матриц, а предлагаете мне прокатиться на чужом не работающем коде. Еще раз пишу вам — ваш вариант не работает, хоть он и быстрее в вашей реализации.


                          1. picul
                            13.10.2019 17:53
                            +1

                            Я дал код для тех, кто пишет в Delphi, а C++ тут даже не упоминается.
                            Если концепция перевода с одного языка на другой Вам недоступна — ОК, вот ссылка на реализацию на Delphi, один в один моя реализация. Я, кстати, возможно даже понял, почему не работает мой код — там сохраняют состояние регистров xmm6-7, чего я не сделал. А вот теперь скажите — даже тем Delphi'стам, кто не может перевести с C++ на Delphi — зачем использовать Ваш медленный вариант, если есть нормальный по выше данной ссылке?
                            Вот и вы пришли совершенно не разбираясь в умножении матриц, а предлагаете мне прокатиться на чужом не работающем коде.
                            То, что Вы не понимаете, как код работает — не означает, что он нерабочий. А перед тем, как делать выводы о том, о чем знает или не знает незнакомый Вам человек — советую почитать про эффект Даннинга-Крюгера, и прочитанным хорошенько проникнуться.
                            Еще раз пишу вам — ваш вариант не работает, хоть он и быстрее в вашей реализации.
                            Еще раз пишу Вам — до Вас в его работоспособности никто не сомневался (а используется он везде, кроме Ваших разработок).


                            1. eDmk Автор
                              13.10.2019 18:07
                              -1

                              Да вы больше на тролля похожи, чем на человека в смысле человек. Мне с вами даже общаться уже нет желания.

                              Код на Delphi по сылке также некорректный, хотя и быстрее вашего. Он повторяет все ошибки вашего некорректного и не рабочего варианта.


                              1. picul
                                13.10.2019 18:18

                                Нда, ну видимо все варианты некорректны, кроме Вашего. Но тролль здесь я, ну так то да, логично.


                                1. truthfinder
                                  14.10.2019 09:16

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


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

                                  Просто восторг.


                      1. mayorovp
                        14.10.2019 08:56
                        +1

                        Я вот не поленился и выполнил все эти действия на бумаге в символьном виде. Матрица транспонировалась корректно.


                        Но писать комментарии про то, что это чушь, конечно же, проще.


                    1. eDmk Автор
                      13.10.2019 17:16

                      Вы мне дайте кол-во умножений вашего общепринятого варианта в секунду. Вам же не составит труда замерить время в цикле? Видимо разговор зашел в тупик.


                      1. picul
                        13.10.2019 17:31

                        Мне не составит труда замерить время в цикле на C++ на моей машине (i5-8250u). Результат — 73 644 864 умножений матриц в секунду.


                        1. eDmk Автор
                          13.10.2019 17:39
                          -1

                          У меня такой же результат +- пару тысяч.
                          Только ваш вариант общепринятый не умножает корректно.


                          1. picul
                            13.10.2019 17:56

                            А в статье написали, что самый быстрый Ваш вариант выдает 52 млн умножений в секунду. И кому теперь верить?
                            И учтите, что Ваш процессор быстрее моего в однопотоке примерно на 20-30%. пруф


                            1. eDmk Автор
                              13.10.2019 18:15

                              Вы из дурдома или тролль? Я вам еще раз говорю вы даете НЕКОРРЕКТНЫЙ КОД и выдаете его за рабочий и общепринятый стандарт. Уйдите отсюда пожалуйста.


                              1. playermet
                                14.10.2019 10:54

                                Этот «некорректный код» используется в том числе в glm, и как следствие в подавляющем большинстве современных графических и игровых движков, приложений и игр на OpenGL/Vulkan и С++. И все как один дают правильную и ожидаемую картинку.


        1. eDmk Автор
          13.10.2019 14:53

          Специально для вас сделал вариант с транспонированной загрузкой.
          Скорость как и говорил немного медленнее: 35 087 719 умножений в секунду.

          Заголовок спойлера
          function MulMatrixCMT(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 из памяти
            //--------------------------------------------
          
            pinsrd xmm4, dword ptr [B + 0], 0
            pinsrd xmm5, dword ptr [B + 4], 0
            pinsrd xmm6, dword ptr [B + 8], 0
            pinsrd xmm7, dword ptr [B + 12], 0
          
            pinsrd xmm4, dword ptr [B + 16], 1
            pinsrd xmm5, dword ptr [B + 20], 1
            pinsrd xmm6, dword ptr [B + 24], 1
            pinsrd xmm7, dword ptr [B + 28], 1
          
            pinsrd xmm4, dword ptr [B + 32], 2
            pinsrd xmm5, dword ptr [B + 36], 2
            pinsrd xmm6, dword ptr [B + 40], 2
            pinsrd xmm7, dword ptr [B + 44], 2
          
            pinsrd xmm4, dword ptr [B + 48], 3
            pinsrd xmm5, dword ptr [B + 52], 3
            pinsrd xmm6, dword ptr [B + 56], 3
            pinsrd xmm7, dword ptr [B + 60], 3
          
            //--------------------------------------------
            // Умножение матриц A?B
            //--------------------------------------------
          
            //------------------------
            // Ряд[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;


          1. picul
            13.10.2019 15:25

            Специально для Вас реализовал нормальный способ на ассемблере:

            Код
            function MulMatrixCMT(A, B: TMatrix4): TMatrix4;
            asm
              .NOFRAME
            
              //--------------------------------------------
              // Загружаем матрицу B строками
              //--------------------------------------------
            
              movups xmm0, dqword ptr [B]       // A[00..03]
              movups xmm1, dqword ptr [B + 16]  // A[04..07]
              movups xmm2, dqword ptr [B + 32]  // A[08..11]
              movups xmm3, dqword ptr [B + 48]  // A[12..15]
              
              //--------------------------------------------
              // Умножение матриц A?B
              //--------------------------------------------
            
              //------------------------
              // Ряд[0]
              //------------------------
              
              movaps xmm4, dqword ptr[A]
              movaps xmm5, xmm4
              movaps xmm6, xmm4
              movaps xmm7, xmm4
            
              shufps xmm4, xmm4, 00000000b
              shufps xmm5, xmm5, 01010101b
              shufps xmm6, xmm6, 10101010b
              shufps xmm7, xmm7, 11111111b
            
              mulps xmm4, xmm0
              mulps xmm5, xmm1
              mulps xmm6, xmm2
              mulps xmm7, xmm3
              
              addps xmm5, xmm4
              addps xmm7, xmm6
              addps xmm7, xmm5
              
              movups dqword ptr [Result], xmm7
              
              //------------------------
              // Ряд[1]
              //------------------------
            
              movaps xmm4, dqword ptr[A + 16]
              movaps xmm5, xmm4
              movaps xmm6, xmm4
              movaps xmm7, xmm4
            
              shufps xmm4, xmm4, 00000000b
              shufps xmm5, xmm5, 01010101b
              shufps xmm6, xmm6, 10101010b
              shufps xmm7, xmm7, 11111111b
            
              mulps xmm4, xmm0
              mulps xmm5, xmm1
              mulps xmm6, xmm2
              mulps xmm7, xmm3
              
              addps xmm5, xmm4
              addps xmm7, xmm6
              addps xmm7, xmm5
              
              movups dqword ptr [Result + 16], xmm7
            
              //------------------------
              // Ряд[2]
              //------------------------
            
              movaps xmm4, dqword ptr[A + 32]
              movaps xmm5, xmm4
              movaps xmm6, xmm4
              movaps xmm7, xmm4
            
              shufps xmm4, xmm4, 00000000b
              shufps xmm5, xmm5, 01010101b
              shufps xmm6, xmm6, 10101010b
              shufps xmm7, xmm7, 11111111b
            
              mulps xmm4, xmm0
              mulps xmm5, xmm1
              mulps xmm6, xmm2
              mulps xmm7, xmm3
              
              addps xmm5, xmm4
              addps xmm7, xmm6
              addps xmm7, xmm5
              
              movups dqword ptr [Result + 32], xmm7
            
              //------------------------
              // Ряд[3]
              //------------------------
            
              movaps xmm4, dqword ptr[A + 48]
              movaps xmm5, xmm4
              movaps xmm6, xmm4
              movaps xmm7, xmm4
            
              shufps xmm4, xmm4, 00000000b
              shufps xmm5, xmm5, 01010101b
              shufps xmm6, xmm6, 10101010b
              shufps xmm7, xmm7, 11111111b
            
              mulps xmm4, xmm0
              mulps xmm5, xmm1
              mulps xmm6, xmm2
              mulps xmm7, xmm3
              
              addps xmm5, xmm4
              addps xmm7, xmm6
              addps xmm7, xmm5
              
              movups dqword ptr [Result + 48], xmm7
              
            end;


            1. eDmk Автор
              13.10.2019 15:39

              Ваш способ быстрее конечно, но он вообще не работает. Будем искать ошибки.


  1. YuraLia
    13.10.2019 13:41

    Была уже отличнейшая статья Умножение матриц: эффективная реализация шаг за шагом на эту же тему, но написана (на мой взгляд) более понятно. Думаю, стоило бы взять с нее пример. Ну и заодно методологию стоило бы повзаимствовать: сравнивать не только в % прирост производительности от оптимизации, но и с теоретической производительностьюю вашего процессора.


    1. eDmk Автор
      13.10.2019 14:25
      -1

      Замечательно, но у меня в теме Delphi + Assembler.
      Вы сравните «эффективную» реализацию по скорости сначала, а потом с советами приходите.
      Я вам свою версию не навязываю.


      1. kenny5660
        13.10.2019 19:13

        А чем Assembler Delphi отличается от assembler'а C++ по вашему? Скорость и возможности ассемблера не зависят от языка высокого уровня который его вызывает. Прочитав всю ветку ваших комментариев, складывается ощущение упертого и высокомерного человека.


        1. eDmk Автор
          13.10.2019 20:49
          -2

          Мне навязывают C++ и нерабочий код и говорят, что это хорошо.
          Где я себя превознес над кем то? Какое тут высокомерие?
          Я же написал, что у меня нет C++, но навязывание не прекратилось.
          Вот и ощущение, что народ или троллит или полностью неадекваты.


          1. playermet
            14.10.2019 10:39

            Вы так много раз назвали метод нерабочим, но может вы просто попробуете взять две реализации и найти хотя бы один набор данных, при которых результат будет отличаться сильней чем на уровне погрешностей? Ну а после десяток/сотен очевидно безуспешных попыток извинитесь перед picul за тон?


            1. eDmk Автор
              14.10.2019 11:18

              Конечно не рабочий. Там нет транспонирования. Уже устал доказывать.
              Извиняться за тон? Я никого не оскорблял, а он упорно навязывает мне свой стандартный общепринятый сишный код. Нормальные люди разбираются в проблеме и находят решение, а не сыплют ссылками с кодом, который мне вообще не нужен.


  1. BD9
    13.10.2019 23:48

    Обращение более к читателям чем к автору:

    1. На Delphi никто расчёты не делает, пользуют C/C++ и Фортран. Можно получить существенный выигрыш просто перейдя на эти языки + соответствующий компилятор, созданный для расчётных задач.
    2. На C/C++ и Фортране можно быстро перейти на GPGPU.
    3. На сегодня без SSE/AVX и прочих SIMD якобы оптимизированный код не нужен.
    4. Используйте библиотечные функции, оптимизированные и отлаженные.
    5. type TVec4 = array[0..3] of Single — использование 32 бит на 64-битной машине — ошибка для расчётчика. Желаете 32 бит и скорость — размещайте 2х32 в одном 64-бит слове.
    6. type TMatrix4 = array[0..3] of TVec4 — для скорости лучше использовать одномерный массив чтобы руками создавать последовательный доступ к памяти. Т.е. двумерный массив преобразуем в одномерный.
    7. type TVec4 = array[0..3] of Single — какой-то это не Паскаль Delphi.


    1. sabaca
      14.10.2019 10:30

      На Delphi никто расчёты не делает, пользуют C/C++ и Фортран.

      Делают, но Вы правы это не лучший выбор.

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

      Ключевое здесь это более качественный оптимизирующий компилятор, который для C/C++/Fortran существуют (и не один), а для Delphi нет (в смысле что качество оптимизации там низкое и нет альтернатив).

      type TVec4 = array[0..3] of Single — использование 32 бит на 64-битной машине — ошибка для расчётчика. Желаете 32 бит и скорость — размещайте 2х32 в одном 64-бит слове.

      Так тут так и получатся.

      type TMatrix4 = array[0..3] of TVec4 — для скорости лучше использовать одномерный массив чтобы руками создавать последовательный доступ к памяти. Т.е. двумерный массив преобразуем в одномерный.

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

      type TVec4 = array[0..3] of Single — какой-то это не Паскаль Delphi.

      Почему?


      1. playermet
        14.10.2019 11:02

        > Почему?
        Думаю, речь о том что в паскалях и дельфи индексы статических массивов обычно с единицы делают.


        1. sabaca
          14.10.2019 11:21

          Ну я всегда с нуля начинал, когда на Паскале писал, там ведь при передаче массива в подпрограмму нижняя граница от нуля начинается.


    1. AndrewBJ
      14.10.2019 11:20

      1. Делают
      2. Возможно.
      3. Что???
      4. Абсолютно верно.
      5. А что не так? Данное описание массивом как раз и размещает в 1 64бит слове 2 по 32. Т.е. элементы в массиве идут непосредственно друг за другом в памяти
      6. Абсолютно всё равно. См п.5.
      7. А что это тогда? Всё корректно написано.


    1. iburyl
      14.10.2019 11:20

      5. type TVec4 = array[0..3] of Single — использование 32 бит на 64-битной машине — ошибка для расчётчика. Желаете 32 бит и скорость — размещайте 2х32 в одном 64-бит слове.

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


  1. AWSVladimir
    14.10.2019 10:21

    1. На Delphi никто расчёты не делает, пользуют C/C++ и Фортран. Можно получить существенный выигрыш просто перейдя на эти языки + соответствующий компилятор, созданный для расчётных задач.

    А можете пример привести?
    Конкретный С/С++ код, конкретный компилятор и его опции в контексте этой статьи?
    Ведь этих Си компиляторов огромная куча, для нужных мне вычислений несколько лет назад в топе был интелловский компилятор, как дело обстоит сейчас?
    PS:
    А на счет Делфи, мне он нравится и так же годится для вычислений, компилятор конечно швах, пальму первенства которая была у Борланда давно уже сломали и сожгли, но использование самому асма или использование в проекте DLL с оптимизированным алгоритмом написанное на любом языке программирования, в плане производительности будет аналогичным.