05.03.2021 12:12
Tiriet 32 3500 Источник

Некоторое время назад понадобилось мне в одной Delphi-шной программе много посчитать, но расчеты шли как-то подозрительно долго. Переписывать около 100 kLOC не хотелось- особенно из-за наличия большого количества форм, а предыдущий мой опыт показывал, что если код расчетов перекомпилировать в Lazarus'е (с FPC3.0.4)- то скорость счета возрастает до 2-х раз, и поэтому было очевидно, что конкретно в данном случае компилятор Embarcadero (разных версий) сильно несилен, и надо его менять. С другой стороны, IDE от Embarcadero для быстрого рисования GUI- вне конкуренции, а их компилятор на редкость быстрый (оно и понятно- быстро+плохо, или медленно+хорошо). Но ведь вкус кактуса неимоверно притягателен. Профайлинг подручными средствами (tinyprofiler) во всех случаях показывал, что основное время (90%) занимают операции векторной алгебры над большими массивами чисел, а быстрый тест производительности этих процедур показал, что на операциях с этой алгеброй общая "пропускная способность" имевшегося математического ядра составляет для операций типа умножения векторов и скалярных произведений- ~4 ГБ/с, для умножения вектора на матрицу- 1,5-2 ГБ/с, а вот для операций обращения матрицы- проседает до 360 МБ/с (на Core I5 4460 и на Xeon 2660V2, DDR3-1866). Внутре рядом с неонкой используются только 3-х и 4-х мерные вектора и матрицы. GUI и ввод-вывод крутятся в своем потоке, а расчеты- крутятся в своих потоках, и друг другу не мешают. В голову пришла мысль, что 4х4 матрица- должна целиком влезать в SSE-регистры процессора и для нее SIMD-очень желательны, а Embarcadero в компилятор SIMD не завезли, не завезут, и вообще- дальше нижней половины XMM0 не используют. В итоге нарисовалась очень простая задача- реализовать быструю векторную алгебру в минимальном объеме для 3D/4D векторов своими руками- то есть, соорудить стероидный велосипед, о котором в заголовке написано.

И вот здесь мой навык гуглинга дал сбой- то есть, навык-то есть- мне предлагали купить проприетарные библиотеки, в которых все уже есть, или примеры векторных операций с FP32 на SSE, но примеров под FP64- не было. Теперь будет. Под катом- как на SSE4.2 руками сделать операцию с вектором для расчетов какой-нибудь физики и вывернуть матрицу наизнанку.

Как я отметил, на тот момент в наличии были Core I5 4460 & Xeon 2660v2, которые поддерживают SSE4.2, но еще не имеют поддержки FMA и AVX, поэтому решать задачу будем на том, что есть, но по ходу пьесы я отмечу те места, где AVX мог бы быть полезен и объясню, почему он бесполезен. Впрочем, объясню сейчас- AVX в контексте бесполезен потому, что даже SSE4.2 позволяет средненькому процу переваривать данные быстрее, чем они пролазят через память, даже одно ядро Xeon'а на 3,3 ГГц обрабатывает ~10ГБ/с, а целый процессор забивает все 4-ре канала памяти и успевает покурить в перерывах, если вся задача не влезает в кэш (у Ryzen 3900 кэша 64МБ, и в такой кеш влезает довольно много, потому на нем смысл в AVX может и возникнуть, но тогда ботлнек случается на переключении между потоками, и заметного выигрыша мне все равно не получить).

Итак, к делу. Минимальные начальные данные:

В нашем процессоре имеется SSE4.2. А это значит, что нам доступны 16 XMM регистров, каждый из которых может содержать по 2 FP64-числа, и всего в XMM мы можем разместить 32 таких числа. Немножко огорчает тот факт, что мы не можем просто так напрямую быстро и эффективно перекидывать данные из XMM в обычные регистры и обратно, но 32 числа- это все-таки довольно много.

Регистры пронумерованы от XMM0 до XMM15. Каждый регистр содержит 128 бит, которые могут рассматриваться по-разному- как просто пачка бит, как 2 FP64, как 4 FP32, как 4 Int32, и много еще как, но нас будет интересовать именно 2 FP64- так как для расчетов физики в моем случае используются только они.

Для таких чисел в SSE4.2 есть ряд очень приятных инструкций, а именно:

movupd XMM4, oWORD[ Pk ] - MOVe Unaligned Packed Double загружает в регистр XMM4 octo-WORD, то есть- 16 байт, из адреса [Pk], причем, загружает их невыровненные, то есть, этот адрес может быть любой, и как бы напоминает нам, что это не просто какие-то байты, а это два вещественных числа.

movapd XMM4, oWORD[ Pk ] - все то же самое, но адрес должен быть выровнен (Aligned) в памяти на 16, если подсунем невыровненный- получим исключение, а если выровненный- то загрузим данные чуть-чуть быстрее, чем movupd.

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

movsd XMM4, qWORD[ Pk ] - MOVe Single Double загружает в нижнюю половину регистра XMM4 quad-WORD, то есть- 8 байт, из адреса [Pk], как бы напоминает- что это именно FP64 число, верхнюю половину не трогает, и к выровненности требований не предъявляет. Можно сделать movsd XMM0, XMM1- тогда она переложит одно нижнее число из XMM1 в XMM0. Ее аналог- movlpd, но она не умеет перекладывать из регистра в регистр, только из/в память.

movlhps XMM0, XMM1 - MOVe Low to High Packed Single - внезапно, берет нижнюю половину XMM1 за два вещественных одинарной точности, и перекладывает их в верхнюю половину XMM0. но так как на самом деле это все просто 64 бита- то получается, что для нас она перекладывает одно двойное число из нижней половины в верхнюю. причем, можно использовать ее в виде movlhps XMM0, XMM0- тогда мы продублируем нижнее число регистра в его же верхнюю половину, и получим два одинаковых числа в одном регистре.

movhlps XMM0, XMM1- наоборот- переложить верхнюю половину одного регистра в нижнюю половину другого (можно и этого же).

mulpd XMM0, XMM4- ооо, это одна из наших любимых инструкций на ближайшее время, MULtiply Packed Doubles- в XMM0 у нас лежат (X0, Y0), в XMM4 у нас лежат (X1, Y1), эта инструкция попарно их перемножает и складывает результат снова в XMM0- то есть, после ее выполнения в XMM0 будет лежать ( X0X1, Y0Y1)- два умножения за время одного!

addpd XMM0, XMM1 - вторая любимая жена инструкция- ADDition Packed Double- после нее XMM0 = ( X0+X1, Y0+Y1).

Иногда бывает еще нужно mulsd/addsd XMM0, XMM1- перемножить/сложить нижние числа регистров и положить результат в нижнюю половину первого операнда.

Еще одна очень полезная инструкция

haddpd XMM0, XMM2 - Horizontal ADDition Packed Double: если XMM0 = (x0,y0), а XMM1 = (x1,y1), то после выполнения этой инструкции XMM0 = (x0+y0, x1+y1)- два сложения за такт по горизонтали (в отличие от addpd, которая складывает "по вертикали").

ну и наверное, последняя из интересных инструкций

dppd XMM0, XMM6,  49 // imm8 = bit0 | ~bit1 | bit4 | bit5 

Dot Product Packed Double - скалярное произведение приемника и источника с маской. Маска выбирает, какие части регистров перемножать, как складывать и куда сохранять результаты. Маска 49- говорит, что нужно X0X1+Y0Y1 положить в нижнюю половину регистра приемника. 

divpd/divsd - как нетрудно догадаться, противоположны mulpd/mulsd и выполняют деление- но они медленные! латентность и время выполнения у них в четыре-пять раз выше.

Так как речь идет о математике, то все функции и процедуры базовых вычислений по умолчанию были с директивами inline и register. Инлайн избавляет от затрат на вызовы в циклах, а register- заставляет компилятор пересылать данные через регистры процессора, а не стек, что несколько быстрее. В документации Delphi сказано, что параметры передаются через регистры, если их меньше пяти, а начиная с пятого- через стек- это надо иметь в виду. Еще в x64-mode Delphi не умеет ассемблерные вставки (только процедуры целиком) и не умеет инлайн ассемблерных функций, так что при переходе на ассемблер придется обходиться без inline- а значит короткие функции будут иметь большой штраф на вызов, и желательно обрабатывать сразу большие пачки данных, чтобы его скомпенсировать :-(.

Хороший лонгрид для начала с ассемблером...

Документация от Embarcadero

Предполагая, что читатель знаком с ассемблером или прочел лонгрид, давайте быстренько перемножим два массива чисел:

T_RealArr = array of real;
T_Mul_s_s = class   
  private    
  fSrc1: T_RealArr;          
  fSrc2: T_RealArr;     
  fRes: T_RealArr;     
  procedure exec_SSE(const i0,i1: integer);   
end;
…

// IMPLEMENTATION
const FloatSize = 8;
procedure T_Mul_s_s.exec_SSE(const i0,i1: integer); 
{$IFDEF CPUX64}
asm   
// Self = RCX   
  .PUSHNV R13
  .PUSHNV R14
  .PUSHNV R15
  XOR R15,R15
  XOR R14,R14
  XOR R13,R13
  mov R13, Self.fRes      // R13 = @Res[0]
  mov R14, Self.fSrc1     // R14 = @Src1[0]
  mov R15, self.fSrc2     // R15 = @Src2[0]

  // расчет смещения Elem[i0] относительно начала массива
  XOR RAX, RAX            // очистили RAX
  mov EAX, i0             // i0- int32, потому кладем ее в EAX, для RAX она маловата
  imul RAX, FloatSize     // а умножать можно весь RAX
  // начальные адреса блоков данных, которые будем обрабатывать.
  add R13, RAX           // R13 = @Res[i0]
  add R14, RAX           // R14 = @Src1[i0]
  add R15, RAX           // R15 = @Src2[i0]

  // а теперь закинем куда-нибудь число оставшихся элементов:
  mov EAX, i1
  sub EAX, i0
  inc EAX // EAX = 1+i1-i0
  
// to be continued…

Массивы у нас лежат в памяти одним куском, но для многопоточности их придется обрабатывать по частям, многопоточностью управляет другой код, а здесь мы запрашиваем диапазон элементов [i0, i1], которые нужно обработать.

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

.PUSHNV R13 говорит компилятору, что перед вызовом нашей функции надо сохранить на стек регистр R13, а после вызова- восстановить его значение, позволяет не париться вопросом- у всякого-ли пуша есть своя попа. для XMM-регистров есть своя, аналогичная директива .SAVENV XMM5

Документация Embarcadero говорит, что при вызове функций регистры XMM0-XMM4 сохраняются и восстанавливаются вызывающим кодом, поэтому об их целостности можно не беспокоиться, а регистры XMM5-XMM15 должны сохраняться и восстанавливаться вызываемым кодом, и поэтому, если мы их используем, то чтобы не поломать что-то снаружи нужно позаботиться об этом, иначе наша процедура может затереть данные, положенные в эти регистры кем-то другим.

mov R13, Self.fRes - требует объяснения- все данные у меня завернуты в объект, Self передается в метод обычно первым параметром через RCX, fRes- имеет какое-то смещение относительно него и его даже можно рассчитать, но Delhpi упрощает мне жизнь, позволяя из своего ассемблера обращаться к переменным и полям по их обычным именам. А так как разные версии компилятора иногда кладут эти переменные по разным регистрам (вопреки документации), то при разработке на нескольких машинах могут возникать странные баги если верить в то, что Self всегда будет в RCХ, а i0- в R8d. Обращение же по имени переменной всегда нормально обрабатывается. Но если придется переносить этот код в FPC- то такое уже не взлетит, придется передавать сами массивы в качестве параметров, да там и вообще синтаксис ассемблера другой.

// continue…
 @new_pair:

    movapd XMM0, oWORD[R14]    // XMM0 = src1[i]   | src1[i+1]
    movapd XMM1, oWORD[R15]    // XMM1 = src2[i]   | src2[i+1]
    mulpd  XMM0, XMM1          // XMM0 = src1[i+0]src2[i+0] | src1[i+1]src2[i+1]
    movapd oWORD[R13], XMM0    // save to Res[i], Res[i+1]

    add R15, FloatSize*2  // next src2 pair
    add R14, FloatSize*2  // next src1 pair
    add R13, FloatSize*2  // next res pair

    sub RAX, 2
    cmp RAX, 1

  JA @new_pair

  cmp RAX, 0

  JE @finish


  @unaligned_end:

      movlpd XMM0, qWORD[R14]    // XMM0 = src1[i]   | garbage0
      movlpd XMM1, qWORD[R15]    // XMM1 = src2[i]   | garbage1
      mulsd  XMM0, XMM1          // XMM0 = src1[i]* src2[i]  | garbage0
      movlpd qWORD[R13], XMM0

      add R15, FloatSize         // move to the next element
      add R14, FloatSize
      add R13, FloatSize

      dec RAX
      cmp RAX,0

  JA @unaligned_end // if RAX > 0- jump to next elem

@finish:

{$ELSE CPUX64} begin {$ENDIF CPUX64}

end;

Cначала мы бежим по нашим данным, загружаем буфер по N=2 числа, перемножаем их тоже пачками, а результаты перекидываем в память и уменьшаем счетчик оставшихся элементов на N. Если останется меньше N элементов- то загрузить полный буфер не удастся, тогда обработаем элементы поштучно, а если останется 0- то мы просто прыгнем на finish. Вместо N=2 легко сделать N=4 (скорости это не добавляет).

На данном коде ускорение по сравнению с обычной (компиляторной) версией примерно в 1.5 раза. Я говорил, что Lazarus мне давал ускорени в 2 раза- но там код брал данные, раскиданные в памяти как попало, и Лазарус почему-то реализовал это сильно эффективнее. Здесь у нас данные лежат в памяти кучно- КМК- компилятору негде особенно косячить, поэтому хуже он просто не смог.

Теперь чуть медленнее попробуем умножить кучу векторов на константу. Так как с выровненными на 16 векторами работать удобнее- то все вектора были превращены в 4-D-вектора, но трехмерный код игнорирует четвертую компоненту (а иногда использует ее для каких-нибудь промежуточных данных, в нее можно скинуть массу частицы или энергию, например, или длину вектора сохранить).

T_Vect = record    
  x: real;     
  y: real;     
  z: real; 
  {$IFDEF ALIGNED_VECTORS_ON}
  t: real;
  {$ENDIF }
end;

T_VectArr = array of T_Vect;
{$IFDEF ALIGNED_VECTORS_ON}
const VectSize = FloatSize * 4;
{$ELSE}
const VectSize = FloatSize * 3;
{$ENDIF }

 T_Mul_v_p = class
  private
    fCoeff: mreal;
    fSrc: T_VectArr;
    procedure exec_SSE(const i0,i1: integer);
  end;

  procedure T_Mul_v_p.exec_SSE_fast(const i0,i1: integer);

{$IFDEF CPUX64}
  asm
  // Self = RCX
  /// Result = Source!
     .NOFRAME
     .PUSHNV R14
     
      XOR RAX, RAX
      mov EAX, i0
      imul RAX, VectSize

      XOR R14, R14
      mov R14, Self.fSrc      // R14 = @Src[0]          
      add R14, RAX            // R14 = @Src[i0]
      
      mov EAX, i1
      sub EAX, i0

      movlpd  XMM2, Self.fCoeff // XMM2 = fCoeff | x
      movlhps XMM2, XMM2        // продублировали в верхнюю половину

  @new_vect:

        movapd XMM0, oWORD[R14]                // XMM0 = v[i].x   | v[i].y
        movapd XMM1, oWORD[R14+2 * FloatSize]  // XMM0 = v[i].z   | v[i].t
        mulpd  XMM0, XMM2
        mulsd  XMM1, XMM2                      // только нижняя половина!

        movapd oWORD[R14              ], XMM0
        {$IFDEF ALIGNED_VECTORS_ON}
        movapd oWORD[R14+2 * FloatSize], XMM1
        {$ELSE}
        movlpd qWORD[R14+2 * FloatSize], XMM1
        {$ENDIF }

        add R14, VectSize  // next vector
        dec RAX
        cmp RAX, 0

  JGE @new_vect // Jump if RAX >=0
{$ENDIF CPUX64}
end;

Код получился почти одинаковый- но это понятно, так как массив векторов- это, фактически, просто в 4-ре раза удлиненный массив скаляров.

Вектора хорошо, но нам нужны еще и тензоры второго ранга в 3-мерном пространстве (тензора деформаций и напряжений, например, или инерции):

T_Tens = packed record
   x: T_Vect;
   y: T_Vect;
   z: T_Vect;
end;

Две часто встречающихся операция с ними- а.i += c.jB.ji и a.i += B.ijc.j- то есть, добавление к вектору правого или левого матричного произведения.

class procedure T_SSE.Add(var a_i: T_Vect; const B_ij: T_Tens; const c_j: T_Vect)
; static; register;

{$IFDEF CPUX64} 
asm 
{   a_i.x := a_i.x+b_ij.x.x * C_j.x+b_ij.x.y * C_j.y+b_ij.x.z * C_j.z ;   
    a_i.y := a_i.y+b_ij.y.x * C_j.x+b_ij.y.y * C_j.y+b_ij.y.z * C_j.z ;   
    a_i.z := a_i.z+b_ij.z.x * C_j.x+b_ij.z.y * C_j.y+b_ij.z.z * C_j.z ;}   
    // C_j: RAX,   
    // b_ij:  RDX, r8   (? compiler dependent! do not use)   
    // a_i: RCX

.NOFRAME
.SAVENV XMM5
.SAVENV XMM6
.SAVENV XMM7

movupd  XMM6, oWORD[C_j+0 * FloatSize]      // XMM6 = c0 | c1
movsd   XMM7, qWORD[C_j+2 * FloatSize]      // XMM7 = c2 | w
movupd  XMM4, oWORD[a_i+0 * FloatSize]      // XMM4 = a0 | a1
movsd   XMM5, qWORD[a_i+2 * FloatSize]      // XMM5 = a2 | w = waste
movupd  XMM0, oWORD[ b_ij+0 * FloatSize ]   //  XMM0 = m00 | m01
movsd   XMM1, qWORD[ b_ij+2 * FloatSize ]   //  XMM1 = m02 | w
movupd  XMM2, oWORD[ b_ij+4 * FloatSize ]   //  XMM2 = m10 | m11
movsd   XMM3, qWORD[ b_ij+6 * FloatSize ]   //  XMM3 = m12 | w
mulpd   XMM0, XMM6            // xmm0 = m00c0 | m01c1
mulsd   XMM1, XMM7            // xmm1 = m02c2 | 0
haddpd  XMM0, XMM1            // xmm0 = m00c0 + m01c1 | m02c2
mulpd   XMM2, XMM6            // xmm2 = m10c0 | m11c1
mulsd   XMM3, XMM7            // xmm3 = m22c2 | 0
haddpd  XMM2, XMM3            // xmm2 = m10c0 + m11c1 | m12c2
haddpd  XMM0, XMM2            // xmm0 = m00c0 + m01c1 + m02c2 | m10c0 + m11c1 + m12c2
addpd   XMM4, XMM0            // xmm4 = a0 + m0jcj | a1 + m1jcj
movupd  XMM2, oWORD[ b_ij+8 * FloatSize  ]   //  XMM2 = m20 | m21
movsd   XMM3, qWORD[ b_ij+10 * FloatSize  ]  //  XMM3 = m22 | w
mulpd   XMM2, XMM6            // XMM2 = m20c0 | m21c1
mulsd   XMM3, XMM7            // XMM3 = m22c2 | 0
movlhps XMM5, XMM3            // XMM5 = a2, m22c2
haddpd  XMM5, XMM2            // xmm3 = m22c2+a2 | m20c0 + m21c1
haddpd  XMM5, XMM5            // xmm5 = m22c2+a2 + m20c0 + m21c1 | m22c2+a2 + m20c0 + m21c1
movupd  oWORD[ a_i ], XMM4
movsd   qWORD[ a_i+2 * FloatSize ], XMM5

{$ELSE} 
begin   
a_i.x := a_i.x+b_ij.x.x * C_j.x+b_ij.x.y * C_j.y+b_ij.x.z * C_j.z;
a_i.y := a_i.y+b_ij.y.x * C_j.x+b_ij.y.y * C_j.y+b_ij.y.z * C_j.z;
a_i.z := a_i.z+b_ij.z.x * C_j.x+b_ij.z.y * C_j.y+b_ij.z.z * C_j.z; 
{$IFEND CPUX64}
end;
class procedure T_SSE.Add(var a_j: T_Vect; const c_i: T_Vect; const B_ij:
    T_Tens); static; register;
 {$IFDEF CPUX64}
asm
{  a_i.x := a_i.x+b_ij.x.xC_j.x+b_ij.x.yC_j.y+b_ij.x.zC_j.z ;
  a_i.y := a_i.y+b_ij.y.xC_j.x+b_ij.y.yC_j.y+b_ij.y.zC_j.z ;
  a_i.z := a_i.z+b_ij.z.xC_j.x+b_ij.z.yC_j.y+b_ij.z.zC_j.z ;}

  // C_i: RDX,
  // b_ij:  r8   (? compiler dependent! do not use)
  // a_j: RCX

    .NOFRAME
    .SAVENV XMM5
    .SAVENV XMM6
    .SAVENV XMM7

    movupd  XMM4, oWORD[ a_j+0 * FloatSize  ]   // XMM4 = a0 | a1
    movsd   XMM5, qWORD[ a_j+2 * FloatSize  ]   // XMM5 = a2 | w
    movapd  XMM0, oWORD[ b_ij  ]     //  XMM0 = m00 | m01
    ADD b_ij, 2 * FloatSize            
    movsd   XMM1, qWORD[ b_ij  ]     //  XMM1 = m02 | w
    ADD b_ij, 2 * FloatSize            
    movsd   XMM7, qWORD[ C_i  ]      //  XMM7 = c0 | 0
    ADD     C_i, FloatSize           
    movlhps XMM7, XMM7               //  XMM7 = c0 | c0
    mulpd   XMM0, XMM7               //  xmm0 = m00c0 | m01c0
    mulsd   XMM1, XMM7               //  xmm1 = m02c0 | 0
    addpd   XMM4, XMM0               //  a0+c0m00 | a1 + c0m01
    addsd   XMM5, XMM1               //  a2+c0m02 | 0
    movapd  XMM0, oWORD[ b_ij  ]     //  XMM0 = m10 | m11
    ADD b_ij, 2 * FloatSize            
    movsd   XMM1, qWORD[ b_ij  ]     //  XMM1 = m12 | w
    ADD b_ij, 2 * FloatSize            
    movsd   XMM7, qWORD[ C_i    ]    // XMM7 = c1 | 0
    ADD     C_i, FloatSize           //  RDX = C_i + 2 elem
    movlhps XMM7, XMM7               // XMM7 = c1 | c1
    mulpd   XMM0, XMM7               // xmm0 = m10c1 | m11c1
    mulsd   XMM1, XMM7               // xmm1 = m12c1 | 0
    addpd   XMM4, XMM0               //  a0+c0m00+m10c1 | a1 + c0m01+m11c1
    addsd   XMM5, XMM1               //  a2+c0m02+m12c1 | 0
    movapd  XMM0, oWORD[ b_ij  ]     //  XMM0 = m20 | m21
    ADD b_ij, 2 * FloatSize            
    movsd   XMM1, qWORD[ b_ij  ]     //  XMM1 = m22 | w
    movsd   XMM7, qWORD[ C_i   ]     // XMM7 = c2 | 0
    movlhps XMM7, XMM7               // XMM7 = c2 | c2
    mulpd   XMM0, XMM7               // xmm0 = m20c2 | m21c2
    mulsd   XMM1, XMM7               // xmm1 = m22c2 | 0
    addpd   XMM4, XMM0               //  a0+c0m00+m10c1+m20c2 | a1 + c0m01+m11c1+m21c2
    addsd   XMM5, XMM1               //  a2+c0m02+m12c1+m22c2 | 0
    movupd  oWORD[a_j   ], XMM4
    movsd   qWORD[a_j+2 * FloatSize], XMM5
{$ELSE}
begin
  a_j.x := a_j.x+b_ij.x.x * C_i.x+b_ij.y.x * C_i.y+b_ij.z.x * C_i.z ;
  a_j.y := a_j.y+b_ij.x.y * C_i.x+b_ij.y.y * C_i.y+b_ij.z.y * C_i.z ;
  a_j.z := a_j.z+b_ij.x.z * C_i.x+b_ij.y.z * C_i.y+b_ij.z.z * C_i.z ;
{$ENDIF CPUX64}

end;

Обратите внимание- во второй процедуре есть код типа ADD C_i, FloatSize, хотя С_i описана как const C_i. Это работает, так как обе процедуры выше описаны как static; register; и имеют меньше 4-х входных параметров, а это значит, что значение C_i передается указателем (потому что record не влазит в регистр), но при этом кладется в регистр, и ее изменение меняет значение указателя в регистре, но не трогает реальные значения в памяти. Такой же фокус можно делать и с константными значениям других типов, которые целиком влезают в регистр- в ассемблерных процедурах они работают как локальные переменные, но нужно учитывать, что если тип переменной целиком помещается в регистр общего назначения- то вместо адреса будет передано само значение, и если тут ошибиться- то дебаг может быть нетривиальным.

Теперь, попробуем вычислить обратную матрицу

Начнем с матриц 3х3. Для вычисления обратной матрицы в этом случае можно использовать два подхода- или использовать прямую формулу- обратная матрица- это транспонированная матрица алгебраических дополнений элементов, деленая на определитель исходной матрицы, или использовать метод Гаусса с расширенной матрицей.

По прямой формуле для матрицы 33

a.xx

a.xy

a.xz

a.yx

a.yy

a.yz

a.zx

a.zy

a.zz

обратная матрица будет равна

+(A.yy*A.zz - A.yz*A.zy)/D

-(A.xy*A.zz - A.xz*A.zy )/D

+(A.xy*A.yz - A.xz*A.yy)/D

-(A.yx*A.zz - A.zx*A.yz)/D

+(A.xx*A.zz - A.xz*A.zx)/D

-(A.xx*A.yz - A.yx*A.xz )/D

+(A.yx*A.zy - A.yy*A.zx)/D

-(A.xx*A.zy - A.xy*A.zx)/D

+(A.xx*A.yy - A.xy*A.yx)/D

D= +(A.yy*A.zz - A.yz*A.zy)*A.xx -(A.xy*A.zz - A.xz*A.zy )*A.yx +(A.xy*A.yz - A.xz*A.yy)*A.zx;

обратим внимание на две вещи- 

  1.  при расчете определителя в скобках стоят те же самые элементы, что и в первой строке конечной матрицы (а множители у этих скобок- это элементы первого столбца исходной матрицы) 

  2. вместо того, чтобы делить на определитель, лучше сначала вычислить 1/D, а потом уже на нее умножать (на Haswell обычная FMUL имеет латентность 5 и требует 1 такт, а FDIV- имеет латентность 10-24 и требует 8-18 тактов! Используемые в SSE mulpd- латентность 5 и 0,5 такта, divsd- латентность 10-20, 8-14 тактов- поэтому везде, где больше одного деления- лучше сначала посчитать 1/D, и потом на нее умножать).

Метод Гаусса с расширенной матрицей предполагает, что мы сначала построим матрицу вида

a.xx

a.xy

a.xz

1

0

0

a.yx

a.yy

a.yz

0

1

0

a.zx

a.zy

a.zz

0

0

1

над этой матрицей проведем преобразования по методу Гаусса, сначала исключая элементы ниже главной диагонали, затем- выше главной диагонали, и после этого получим матрицу вида:

1

0

0

b.xx

b.xy

b.xz

0

1

0

b.yx

b.yy

b.yz

0

0

1

b.zx

b.zy

b.zz

полученная в правой части b- это и есть обратная матрица.

Экспериментально проверено: для матриц 3х3 прямой расчет через дополнения выполняется быстрее исключения Гауссом. Кроме того, в методе Гаусса могут возникать нули на главной диагонали, что требует переупорядочивания строк, большого количества условных переходов, проверок, и еще усложняет без того не очень эффективный алгоритм, а прямая формула использует 1 деление, требует только одной проверки (Abs(D)> eps), и потому предпочтительнее.

Отметим еще один момент: сама матрица содержит 9 элементов, матрица результатов- еще 9 элементов, а нам доступно 32 элемента, то есть, все расчеты полностью влезают в регистры и остается еще свободное место на хранения промежуточных значений. 

Страшный ассемблерный код (с комментариями), реализующий прямой расчет обратной матрицы 3х3 приведен ниже:

class procedure T_SSE.Invert(const S: T_Tens; const i0, i1: integer);
{$IFDEF CPUX64}
  const One: double = 1.0;
        Zero: double = 0.0;
        eps: double = 1.0e-10;
        sign: Uint64 = $7FFFFFFFFFFFFFFF;
  asm
    .NOFRAME
    .SAVENV XMM5   // S.z.z
    .SAVENV XMM6   // Res.x.x | Res.x.y
    .SAVENV XMM7   // Res.x.z | w
    .SAVENV XMM8   // Res.y.x | Res.y.y
    .SAVENV XMM9   // Res.y.z | w
    .SAVENV XMM10  // Res.z.x | Res.z.y
    .SAVENV XMM11  // Res.z.z | w
    .SAVENV XMM12  // tmp1
    .SAVENV XMM13  // tmp2        
    XOR RAX, RAX
    mov EAX, i0
    imul RAX, 3*VectorSize // Size of Tensor
    ADD RCX, RAX // RCX = @S[i0]
    XOR RAX, RAX
    mov EAX, i1
    sub EAX, i0
    inc RAX
    PXOR XMM6, XMM6
    PXOR XMM7, XMM7
    PXOR XMM8, XMM8
    PXOR XMM9, XMM9
    PXOR XMM10, XMM10
    PXOR XMM11, XMM11
@begin:
    prefetch [RCX+ 12*FloatSize  ]	   // предзагрузка следующего тензора в кэш
    prefetch [RCX+ 16*FloatSize  ]      
    prefetch [RCX+ 20*FloatSize  ]
    
    movupd XMM0, [RCX              ]   //   S.x.x | S.x.y
    movupd XMM1, [RCX+  2*FloatSize]   //   S.x.z | w
    movupd XMM2, [RCX+  4*FloatSize]   //   S.y.x | S.y.y
    movupd XMM3, [RCX+  6*FloatSize]   //   S.y.z | w
    movupd XMM4, [RCX+  8*FloatSize]   //   S.z.x | S.z.y
    movupd XMM5, [RCX+ 10*FloatSize]   //   S.z.z | w
  // Рассчет алгебраических дополнений (с учетом знака перестановки)  
  {   Res.x.x :=    ( + S.y.y*S.z.z - S.y.z*S.z.y  );   
      Res.x.y :=    ( + S.x.z*S.z.y - S.x.y*S.z.z  );          }
    movhlps XMM6, XMM2     // XMM6.L = S.y.y
    movlhps XMM6, XMM3     // XMM6.h = S.y.z
    blendpd XMM5, XMM4, 2  // XMM5= S.z.z | S.z.y
    mulpd XMM6, XMM5       // XMM6 = S.y.yS.z.z |  S.y.zS.z.y
    movhlps XMM12, XMM4    // XMM12.L = S.z.y
    movlhps XMM12, XMM5    // XMM12.h = S.z.z
    blendpd XMM1,  XMM0, 2 // XMM1 = S.x.z | S.x.y
    mulpd XMM12, XMM1      // XMM13 = S.x.z*S.z.y | S.x.y*S.z.z
    hsubpd XMM6, XMM12     // XMM6 = Res.x.x | Res.x.y     }
  {  Res.y.x :=    ( + S.z.x*S.y.z - S.y.x*S.z.z  );
     Res.y.y :=    ( + S.x.x*S.z.z - S.x.z*S.z.x  );   }
    movsd   XMM8, XMM4     // XMM8.L = S.z.x
    movlhps XMM8, XMM2     // XMM8.h = S.y.x
    movlhps XMM3 , XMM5    // XMM3 = S.y.z | S.z.z
    mulpd   XMM8, XMM3     // XMM8 = S.z.x*S.y.z |  S.y.x*S.z.z
    movsd   XMM13, XMM0    // XMM12.L = S.x.x
    movlhps XMM13, XMM1    // XMM13 = S.x.x | S.x.z
    movlhps XMM5, XMM4     // XMM5 = S.z.z  | S.z.x
    mulpd XMM13, XMM5      // XMM13 = S.x.x*S.z.z | S.x.z*S.z.x
    hsubpd XMM8, XMM13     // XMM8 = Res.y.x | Res.y.y
  {   Res.z.x :=    ( + S.y.x*S.z.y - S.y.y*S.z.x  );
      Res.z.y :=    ( + S.x.y*S.z.x - S.x.x*S.z.y  );   }
    movhlps XMM10, XMM4
    movlhps XMM10, XMM4    // XMM10 = S.z.y | S.z.x
    pxor XMM12, XMM12
    subpd  XMM12, XMM10    // XMM10 = -S.z.y | -S.z.x
    mulpd  XMM12, XMM0     // XMM12 = S.z.y*S.y.x | S.z.x*S.y.y
    mulpd  XMM10, XMM2     // XMM10 = -S.x.x*S.z.y | -S.x.y*S.z.x
    hsubpd XMM10, XMM12    // XMM10 = Res.z.x | Res.z.y
  {  Res.x.z :=    ( + S.x.y*S.y.z - S.x.z*S.y.y  ); //  XMM7.L
     Res.y.z :=    ( + S.y.x*S.x.z - S.x.x*S.y.z  ); // XMM9.L
    // Reorder:
    Res.x.z :=    ( + S.x.y*S.y.z - S.y.y*S.x.z  ); //  XMM7.L
    Res.y.z :=   -(   S.x.x*S.y.z - S.y.x*S.x.z  ); // XMM9.L       }
    movlhps XMM3, XMM3     // XMM3 = S.y.z | S.y.z
    mulpd   XMM3, XMM0     // XMM3 = S.x.x*S.y.z | S.x.y*S.y.z
    movapd  XMM12, XMM2    // XMM12 = S.y.x | S.y.y
    movlhps XMM1, XMM1     // XMM1  = S.x.z | S.x.z
    mulpd   XMM12, XMM1    // XMM12 = S.y.x*S.x.z | S.y.y*S.x.z
    subpd   XMM3, XMM12    // XMM3 = S.x.x*S.y.z-S.y.x*S.x.z | S.x.y*S.y.z - S.y.y*S.x.z
                           //      = -Res.y.z                | Res.x.z
    pxor  XMM9, XMM9       // XMM9 = 0 | 0
    subsd XMM9, XMM3       // XMM9 = Res.y.z
    movhlps XMM7, XMM3     // XMM11 = Res.x.z
  {  Res.z.z :=    ( + S.x.x*S.y.y - S.x.y*S.y.x  );     }
    movhlps XMM11, XMM2    // XMM10 = S.y.y | w
    movlhps XMM11, XMM2    // XMM10 = S.y.y | S.y.x
    mulpd   XMM11, XMM0    // XMM11 = S.x.x*S.y.y | S.x.y*S.y.x
    hsubpd  XMM11, XMM11   // XMM11 = S.x.x*S.y.y - S.x.y*S.y.x
  {   D := S.x.x*Res.x.x +   S.x.y*Res.y.x +   S.x.z*Res.z.x ==
           Res.x.x*S.x.x +   Res.x.y*S.y.x +   Res.x.z*S.z.x         }
    movlhps XMM0, XMM2
    dppd    XMM0, XMM6,  49 // dppd- dot product packed double
    												// imm8 = bit0 | ~bit1 | bit4 | bit5
                            // bit0- save res to low half
                            // bit1- save res to high half
                            // bit4- * low parts
                            // bit5- * high parts
                            // XXM0 = (s.xx | s.yx)*(Res.xx | Res.xy)
    mulsd XMM4, XMM7        // XMM4 =  S.zx * Res.xz
    addsd XMM4, XMM0        // now XMM4 = D = Determinant(S)!
    // Check if D< eps:	
    movsd XMM0, XMM4        
    movsd XMM1, sign
    pand  XMM0, XMM1        // убрали знак D- получили XMM0 = Abs(D)
    movsd XMM1, eps
    comisd XMM1, XMM0       // if ABS (D) > eps then PROCEED else EXIT
    jnbe @exit
  //  Res now is the Matrix of algebraic complements, 
  //  Out := mult(Res,1.0/D);
    movsd   XMM3, One
    divsd   XMM3, XMM4
    movlhps XMM3, XMM3	                     // XMM3 = 1.0/D | 1.0/D
    mulpd   XMM6, XMM3
    movupd  oWORD[RCX              ], XMM6   //   S.x.x | S.x.y
    mulpd   XMM7, XMM3
    movupd  oWORD[RCX+  2*FloatSize], XMM7   //   S.x.z |
    mulpd   XMM8, XMM3
    movupd  oWORD[RCX+  4*FloatSize], XMM8   //   S.y.x | S.y.y
    mulpd   XMM9, XMM3
    movupd  oWORD[RCX+  6*FloatSize], XMM9   //   S.y.z |
    mulpd   XMM10, XMM3
    movupd  oWORD[RCX+  8*FloatSize], XMM10  //   S.z.x | S.z.y
    mulpd   XMM11, XMM3
    movupd  oWORD[RCX+ 10*FloatSize], XMM11  //   S.z.z |
  //  sfence
    @exit:
    Add RCX, 3*VectorSize  // next matrix
    dec RAX
    cmp RAX, 0
    jg @begin              // if RAX > 0- invert next matrix
{$ELSE}
var i: integer;
begin
  for i := i0 to i1 do
    Invert(  T_TensorArr(@S)[i] );
{$IFEND}
end;

Хорошо, обращение матрицы 3х3- вещь полезная и часто нужная, но иногда еще нужно обращать матрицы 4х4 (например, при расчете линейной аппроксимации трехмерной функции методом наименьших квадратов на неструктурированном множестве точек). 

Матрицы 4х4

И тут ситуация в корне изменилась.

Во-первых, прямой метод требует расчета алгебраических дополнений- сейчас это стали определители 3х3.

Во-вторых, сама матрица занимает уже 16 элементов- ровно половину доступных нам регистров, матрица результата- еще 16 элементов, промежуточные значения хранить уже физически негде, а перекидывать данные из XMM-регистров в обычные и обратно мы не можем. Можно скидывать данные в стек, но тогда реализация станет совсем неподьемной. Что же делать?

К счастью, у меня созрел гениальный план. Во-первых, в моей физике используемые матрицы- обычно достаточно хорошо обусловлены и как правило- имеют диагональное если не преобладание, то как минимум- превалирование (преобладание- это когда элемент больше суммы модулей остальных элементов, превалирование- это когда он больше по модулю любого другого в строке, но не их суммы). Во-вторых- я подумал, что мне совсем не обязательно хранить всю расширенную матрицу, если я использую метод Гаусса.

Дело в том, что в методе Гаусса с расширенной матрицей на любом этапе очень много фиксированных значений. Вот смотрите- до начала расчетов :

[xx

xy]

[xz

xt]

[1

0]

[0

0]

[yx

yy]

[yz

yt]

[0

1]

[0

0]

[zx

zy]

[zz

zt]

[0

0]

[1

0]

[tx

ty]

[tz

tt]

[0

0]

[0

1]

фиксированные значения я выделил жирным. После зануления первого столбца:

[1

xy’]

[xz’

xt’]

[res.xx

0

[0

0]

[0

yy’]

[yz’

yt’]

[res.yx

1

[0

0]

[0

zy’]

[zz’

zt’]

[res.zx

0

[1

0]

[0

ty’]

[tz’

tt’]

[res.tx

0

[0

1]

с точки зрения операций над строками этой матрицы, нам ничто не мешает перекинуть 5-й столбик вместо первого, но тогда у нас снова окажется в левой половине матрица с какими-то значениями, а в правой- та же самая единичная матрица, с которой мы начали- как будто ничего с ней и не произошло. А раз так- то зачем же ее обрабатывать? она же не меняется!

после зануления элементов ниже главной диагонали мы получим матрицу:

1

xy

xz

xt

res.xx’

0

0

0

0

1

yz

yt

res.yx

res.yy

0

0

0

1

1

zt

res.zx

res.zy

res.zz

0

0

0

0

1

res.tx

res.ty

res.tz

res.tt

И на любом шаге алгоритма ровно половина элементов имеют точно известные значения (нули и единицы). А это значит, нам на самом деле не нужно хранения двух полных матриц, достаточно хранить только изменяющиеся значения. Кроме того, можно заметить, что элементы в левой и правой частях расширенной матрицы очень удачно дополняют друг друга- если слева- полная матрица- то справа- нет вообще ни одного произвольного элемента, если слева- нижний треугольник занулен- то в правой части именно в этом треугольнике и есть ненулевые элементы. Если на главной диагонали левой части стоят единицы- то на главной диагонали правой части- стоят значения, а если в левой части- значения, то в правой части- единицы. Поэтому можно крутить данные только в одной квадратной матрице, используя для хранения промежуточных вычислений 8 освободившихся XMM-регистров. 

Итак, алгоритм у нас получается примерно такой: 

Выбираем строку с номером i=1.
сохраняем D.ii с главной диагонали этой строки, например- “xx”, 
на его место записываем 1.0. 
вычисляем 1/D, 
всю строку с этим элементом умножаем на 1/D. 
	выбираем j-ю строку j=i+1
	выбираем в ней i-й элемент (‘yx’), копируем его.
	вместо него записываем 0. 
	вычисляем 1/’yx’
	умножаем на это значение всю j-ю строку, 
	вычитаем из строки j строку i. 
	переходим к j+1
	если строки кончились- увеличиваем i=i+1, и повторяем. 
когда обработаем так все строки- аналогичным образом выполняем обратный ход метода Гаусса

Из вспомогательных значений нам нужны 0.0e+0, 1.0e+0 (их можно хранить в 1 XMM-регистре), какой-то регистр для копирования промежуточных значений, регистр, в который мы будем сохранять результаты деления, и у нас остается еще 5 свободных регистров, которые можно использовать (но некуда).

Соответствующая портянка ассемблерного кода- под катом. 

type
  T_M4 = array [0..3,0..3] of real;
class procedure T_SSE.Invert_gauss(var S: T_M4);
{$IFDEF CPUX64}
const FloatSize = 8;
      One: double = 1.0;
      Zero: double = 0.0;
asm
  .NOFRAME
  .SAVENV XMM5
  .SAVENV XMM6
  .SAVENV XMM7
  .SAVENV XMM8
  .SAVENV XMM9
  .SAVENV XMM10
  movupd XMM0, [RCX              ]   //   S.x.x | S.x.y
  movupd XMM1, [RCX+  2*FloatSize]   //   S.x.z | S.x.t
  prefetch [RCX + 4*FloatSize]       // slightly increase overall performance (~1%)
  PXOR XMM10,XMM10
//  FORWARD MOVE 
 // STEP  #1
 // STEP #1.prep:  D = 1.0/ S.x.x    S.x.x = 1; S.x := S.xD;
  movsd XMM8, One             // XMM8.L = 1.0;
  movsd XMM9, XMM0            // EXTRACT DIVIDER!
  movsd XMM0, XMM8            //XMM0 = 1 | S.x.y
  divsd XMM8, XMM9            // D := 1.0 / S.x.x;
  movlhps XMM8, XMM8          // XMM8 = D  | D
  mulpd XMM0, XMM8            // S.x = S.xD
  mulpd XMM1, XMM8            // XMM0 XMM1= a00d | a01d   ||  a02d | a03d
  movupd XMM2, [RCX+  4*FloatSize]   //   S.y.x | S.y.y
  movupd XMM3, [RCX+  6*FloatSize]   //   S.y.z | S.y.t
  prefetch [RCX + 8*FloatSize]
 // STEP  #1.y:   D := S.y.x;  s.y.x :=  0; S.y = S.y- S.xD
  movsd   XMM8, XMM2   // XMM8 = (D= S.y.x) | 0
  movlhps XMM8, XMM8   // XMM8 =  D          |   D
  movsd XMM2, XMM10    // S.y.x := 0;  XMM2 = 0 | S.y.y
  movapd XMM9, XMM0    // copy S.x.x, S.x.y
  mulpd XMM9, XMM8     // mul
  subpd XMM2, XMM9     // XMM2 = -S.x.xD         | S.y.y - S.x.yD
  movapd XMM9, XMM1    // copy S.x.z, S.x.t
  mulpd XMM9, XMM8     // mul
  subpd XMM3, XMM9     // XMM3 = S.y.z - S.x.zd |  S.y.t - S.x.td        //}
 // STEP  #1.z:   D := S.z.x;  s.z.x :=  0; S.z = S.z- S.xD
  movupd XMM4, [RCX+  8*FloatSize]   //   S.z.x | S.z.y
  movupd XMM5, [RCX+ 10*FloatSize]   //   S.z.z | S.z.t
  prefetch [RCX + 12*FloatSize]
  movsd   XMM8, XMM4   // XMM8 = (D= S.y.x) | 0
  movlhps XMM8, XMM8   // XMM8 =  D          |   D
  movsd XMM4, XMM10    // S.z.x:= 0;  XMM4 = 0 | S.z.y
  movapd XMM9, XMM0    // Copy S.x.x, S.x.y
  mulpd XMM9, XMM8     // mul
  subpd XMM4, XMM9     // XMM4 = -S.x.xD         | S.y.y - S.x.yD
  movapd XMM9, XMM1    // copy S.x.z , S.x.t
  mulpd XMM9, XMM8     // mul
  subpd XMM5, XMM9     // XMM5 = S.z.z - S.z.zd |  S.z.t - S.z.td        //}
 // STEP  #1.t:   D := S.t.x;  s.t.x :=  0; S.t = S.t- S.xD
  movupd XMM6, [RCX+ 12*FloatSize]   //   S.t.x | S.t.y
  movupd XMM7, [RCX+ 14*FloatSize]   //   S.t.z | S.t.t
  movsd   XMM8, XMM6   // XMM8 = (D= S.y.x) | 0
  movlhps XMM8, XMM8   // XMM8 =  D          |   D
  movsd XMM6, XMM10    // S.t.x:= 0;  XMM6 = 0 | S.z.y
  movapd XMM9, XMM0    // Copy S.y.x, S.y.y
  mulpd XMM9, XMM8     // mul
  subpd XMM6, XMM9     // XMM6 = -S.x.xD         | S.y.y - S.x.yD
  movapd XMM9, XMM1    // copy S.z.z , S.z.t
  mulpd XMM9, XMM8     // mul
  subpd XMM7, XMM9     // XMM7 = S.z.z - S.z.zd |  S.z.t - S.z.td        //}
// STEP  #4.prep  D := 1.0/S.y.y;    S.y.y = 1;  S.y := S.yD;
  movsd   XMM8, One
  movhlps XMM9, XMM2    // XMM9 = S.y.y    |    w
  movlhps XMM2, XMM8    // XMM2 = S.y.x    |    S.y.y=1
  divsd   XMM8, XMM9    // XMM8 = 1.0/ S.y.y   = D
  movlhps XMM8, XMM8
  mulpd   XMM2, XMM8    // XMM2 = S.y.xD  |    D= 1/S.y.y
  mulpd   XMM3, XMM8    // XMM3 = S.y.zD  |    S.y.tD
// STEP  #4.z  D := S.z.y; S.z.y := 0; S.z := S.z-S.yD
  movhlps XMM8, XMM4   // XMM9= D=S.z.y   |   w
  movlhps XMM8, XMM8   // XMM8= D         |   D
  movlhps XMM4, XMM10  // XMM4= S.z.x     | S.z.y=0
  movapd  XMM9, XMM2   // XMM9= S.y.x      | S.y.y
  mulpd   XMM9, XMM8   // XMM9= S.y.xD    | S.y.yD
  subpd   XMM4, XMM9   // XMM4= S.z.x - S.y.xD | 0 -S.y.yD
  movapd  XMM9, XMM3   // XMM9= S.y.z      | S.y.t
  mulpd   XMM9, XMM8   // XMM9= S.y.zD    | S.y.tD
  subpd   XMM5, XMM9   // XMM4= S.z.z - S.y.zD | S.z.t - S.y.tD
// STEP  #4.t  D := S.t.y; S.t.y := 0; S.t := S.t-S.yD
  movhlps XMM8, XMM6   // XMM9= D=S.t.y   |   w
  movlhps XMM8, XMM8   // XMM8= D         |   D
  movlhps XMM6, XMM10  // XMM4= S.t.x     | S.t.y=0
  movapd  XMM9, XMM2   // XMM9= S.y.x      | S.y.y
  mulpd   XMM9, XMM8   // XMM9= S.y.xD    | S.y.yD
  subpd   XMM6, XMM9   // XMM4= S.t.x - S.y.xD | 0 -S.y.yD
  movapd  XMM9, XMM3   // XMM9= S.y.z      | S.y.t
  mulpd   XMM9, XMM8   // XMM9= S.y.zD    | S.y.tD
  subpd   XMM7, XMM9   // XMM4= S.t.z - S.y.zD | S.t.t - S.y.tD
// STEP  #5.prep  D := 1.0/S.z.z;    S.z.z = 1;  S.z := S.zD;
  movsd   XMM8, One
  movsd   XMM9, XMM5    // XMM9 = S.z.z    |    w
  movsd   XMM5, XMM8    // XMM5 = S.z.z=1  |    S.z.t
  divsd   XMM8, XMM9    // XMM8 = 1.0/ S.z.z   = D
  movlhps XMM8, XMM8
  mulpd   XMM4, XMM8    // XMM4 = S.z.xD  |    D= 1/S.z.z
  mulpd   XMM5, XMM8    // XMM5 = S.z.zD  |    S.z.tD
// STEP  #5.t  D := S.t.z; S.t.z := 0; S.t := S.t-S.zD
  movsd   XMM8, XMM7   // XMM8= D=S.t.z   |   w
  movlhps XMM8, XMM8   // XMM8= D         |   D
  movsd   XMM7, XMM10  // XMM7= S.t.z=0   | S.t.t
  movapd  XMM9, XMM4   // XMM9= S.z.x      | S.z.y
  mulpd   XMM9, XMM8   // XMM9= S.z.xD    | S.z.yD
  subpd   XMM6, XMM9   // XMM4= S.t.x - S.z.xD | S.t.y -S.z.yD
  movapd  XMM9, XMM5   // XMM9= S.z.z      | S.z.t
  mulpd   XMM9, XMM8   // XMM9= S.z.zD    | S.z.tD
  subpd   XMM7, XMM9   // XMM7= 0 - S.z.zD | S.t.t - S.z.tD
// STEP  #6  D := 1/S.t.t; S.t.t := 1.0; S.t := S.tD
  movsd   XMM8, One
  movhlps XMM9, XMM7    // XMM9 = S.t.t    |    w
  movlhps XMM7, XMM8    // XMM7 = S.t.z    |    S.t.t=1
  divsd   XMM8, XMM9    // XMM8 = 1.0/ S.z.z   = D
  movlhps XMM8, XMM8
  mulpd XMM6, XMM8
  mulpd XMM7, XMM8
  // save S.t
  movupd oWORD[RCX+ 12*FloatSize], XMM6   //   S.z.x | S.z.y
  movupd oWORD[RCX+ 14*FloatSize], XMM7   //   S.z.z |
// BACKWARD MOVEMENT!
// STEP  #7.x-y  D := S.x.y; S.x.y := 0; S.x -= S.yD
  movhlps XMM8, XMM0     // XMM8= D=S.x.y | w
  movlhps XMM8, XMM8     // XMM8 = D | D
  movlhps XMM0, XMM10    // S.x.y = 0
  movapd XMM9, XMM2      // XMM9= S.y.x      | S.y.y
  mulpd XMM9, XMM8       // XMM9= S.y.xD    | S.y.yD
  subpd XMM0, XMM9       // XMM0= S.x.x+S.y.xD | S.y.yD
  movapd XMM9, XMM3       // XMM9= S.y.z      | t
  mulpd XMM9, XMM8       // XMM9= S.y.zD    | t
  subpd XMM1, XMM9       // XMM0= S.x.z+S.y.zD | t
// STEP  #7.y-z  D := S.y.z; S.y.z := 0; S.y -= S.zD
  movsd   XMM8, XMM3     // XMM8= D=S.y.z | w
  movsd   XMM3, XMM10    // S.y.z := 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM4      // XMM9= S.z.x      | S.z.y
  mulpd XMM9, XMM8       // XMM9= S.z.xD    | S.z.yD
  subpd XMM2, XMM9       // XMM0= S.y.x-S.z.xD | S.y.y-S.z.yD
  movapd XMM9, XMM5       // XMM9= S.y.z      | t
  mulpd XMM9, XMM8       // XMM9= S.y.zD    | t
  subpd XMM3, XMM9       // XMM0= S.x.z+S.y.zD | t
// STEP  #7.x-z  D := S.x.z; S.x.z := 0; S.x -= S.zD
  movsd   XMM8, XMM1     // XMM8= D=S.x.z | w
  movsd   XMM1, XMM10    // S.x.z := 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM4      // XMM9= S.t.x      | S.t.y
  mulpd XMM9, XMM8       // XMM9= S.t.xD    | S.t.yD
  subpd XMM0, XMM9       // XMM0= S.z.x-S.t.xD | S.z.y-S.t.yD
  movapd XMM9, XMM5      // XMM9= S.t.z      | S.t.t
  mulpd XMM9, XMM8       // XMM9= S.t.zD    | S.t.tD
  subpd XMM1, XMM9       // XMM0= S.z.z-S.t.zD | S.z.t-S.t.tD
// STEP  #7.z-t  D := S.z.t; S.z.t := 0; S.z -= S.tD
  movhlps XMM8, XMM5     // XMM8= D=S.z.t | w
  movlhps XMM5, XMM10    // S.z.t := 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM6      // XMM9= S.t.x      | S.t.y
  mulpd XMM9, XMM8       // XMM9= S.t.xD    | S.t.yD
  subpd XMM4, XMM9       // XMM0= S.z.x-S.t.xD | S.z.y-S.t.yD
  movapd XMM9, XMM7      // XMM9= S.t.z      | S.t.t
  mulpd XMM9, XMM8       // XMM9= S.t.zD    | S.t.tD
  subpd XMM5, XMM9       // XMM0= S.z.z-S.t.zD | S.z.t-S.t.tD
  // SAVE S.z
  movupd oWORD[RCX+  8*FloatSize], XMM4   //   S.z.x | S.z.y
  movupd oWORD[RCX+ 10*FloatSize], XMM5   //   S.z.z |
// STEP  #8.y-t  D := S.y.t; S.y.t := 0; S.y -= S.tD
  movhlps XMM8, XMM3     // XMM8= D=S.y.t | w
  movlhps XMM3, XMM10    // S.y.t = 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM6      // XMM9= S.y.x      | S.y.y
  mulpd XMM9, XMM8       // XMM9= S.y.xD    | S.y.yD
  subpd XMM2, XMM9       // XMM0= S.x.x+S.y.xD | S.y.yD
  movapd XMM9, XMM7      // XMM9= S.y.z      | t
  mulpd XMM9, XMM8       // XMM9= S.y.zD    | t
  subpd XMM3, XMM9       // XMM0= S.x.z+S.y.zD | t
//   save S.y
  movupd oWORD[RCX+  4*FloatSize], XMM2   //   S.y.x | S.y.y
  movupd oWORD[RCX+  6*FloatSize], XMM3   //   S.y.z |
// STEP  #8.x-t  D := S.x.t; S.x.t := 0; S.x -= S.tD
  movhlps XMM8, XMM1     // XMM8= D=S.x.t | w
  movlhps XMM1, XMM10    // S.x.t := 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM6      // XMM9= S.t.x      | S.t.y
  mulpd XMM9, XMM8       // XMM9= S.t.xD    | S.t.yD
  subpd XMM0, XMM9       // XMM0= S.x.x-S.t.xD | S.x.y-S.t.yD
  movapd XMM9, XMM7      // XMM9= S.t.z         | S.t.t
  mulpd XMM9, XMM8       // XMM9= S.t.zD       | S.t.tD
  subpd XMM1, XMM9       // XMM0= S.x.z-S.t.zD | S.x.t-S.t.tD
  movupd oWORD[RCX              ], XMM0   //   S.x.x | S.x.y
  movupd oWORD[RCX+  2*FloatSize], XMM1   //   S.x.z |
 {$ELSE}
begin
  Invert_m4( S );
 {$ENDIF}
end;

Здесь по коду раскиданы инструкции вида

 prefetch [RCX + 4FloatSize] 

данная команда предлагает процессору сделать упреждающую выборку данных из памяти, выполняется параллельно с последующим потоком команд, и при возможности- загружает в кэш процессора данные, которые понадобятся чуть позже, в результате доступ к ним будет происходить чуточку быстрее, но максимум, что я заметил с этих префетчей- это 1-2% ускорения. Еще здесь много операций, которые можно было бы заменить на FMA, которая есть в Ryzen и новых интелях, но после проведения ряда тестов желание использовать FMA у меня пропало.

Тесты

На момент написания у меня в наличии уже Xeon-2678v3 c DDR3-1866-4x (без анлока турбобуста- почему- в самом конце) и Ryzen 3900x c DDR4-3200-2x. Тесты будем проводить на Xeon-е, а часть интересных особенностей сравним с Ryzen-ом.  Скорость обработки данных определим как объем данных, который мы можем считать и обсчитать за 1 секунду (без учета того, сколько при этом будет записано результатов в память). Все тесты запускались по 5 раз и выбиралось наилучшее значение. Во всех случаях использовались массивы по 1 048 576 элементов.

Сумма скалярных произведений элементов двух массивов случайных векторов:

Dot product: compiler :  5240,8 MB/s
Dot product: SSE      :  7717,2 MB/s
Dot product: SSE, arr :  9704,8 MB/s

Рост производительности просто от использования SSE небольшой, потому что много времени занимает вызов функции SSE-умножения, так как Embarcadero не умеет инлайн ассемблерных функций.  “Dot product: SSE, arr”- это то же самое умножение у которого цикл тоже реализован на ассемблере внутри SSE-блока, как видим, избавившись от вызовов процедуры мы сократили время расчетов на четверть! Нужно учесть, что умножаются два массива векторов, а результат сохраняется в массив вещественных чисел. В данном бенчмарке мы считываем два раза по 4.8 GB/s данных, и потом еще 1200 MB/s записываем, в сумме прокачивая через канал памяти ~11 GB/s с одного ядра. 

Умножение вектора на матрицу с суммированием результатов:

V3 += V3M33 compiler : 3599,5 MB/s
V3 += V3M33 SSE      : 6810,8 MB/s
V3 += M33V3 compiler : 4304,6 MB/s
V3 += M33V3 SSE      : 8644,2 MB/s

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

M3х3: Invert compiler    : 3384,8 MB/s
M3х3: T_SSE.Invert_gauss : 3386,8 MB/s
M3х3: T_SSE.Invert       : 4990,0 MB/s

Как видим, в моей реализации метод Гаусса на SSE также хорош (или так же плох), как и прямой расчет от компилятора, но вот прямой расчет на SSE- почти в полтора раза быстрее. 

M4х4: Invert compiler    :  333,2 MB/s
M4х4: Invert_SSE_gauss   :  2958,2 MB/s

А вот здесь уже разница стала подавляющей- ускорение почти в 9 раз! Так как исходная программа довольно много аппроксимировала значения МНК- то для нее расчет обратной матрицы 4х4 был самой тормозной операцией, которая занимала почти 60% времени. после оптимизации она стала занимать порядка 10%, а общая скорость счета выросла примерно в 3 раза. 

С матрицам интересно еще вот что: зависимость количества операций, необходимых для инверсии матрицы, от ее размера имеет характер O(n5), то есть, при переходе с М3x3 на M4x4- объем вычислений возрастает в ~4,2 раза, объем данных при этом возрастает в 1,8 раза, скорость обработки данных должна упасть в 2,3 раза, однако, падает только в 1.7 раза, то есть, меньше, чем следует. А связано это с тем, что данные в регистрах лежат плотнее (у M3x3 строка занимает полтора регистра, а у M4x4- ровно два) , и параллельность их обработки в целом получается выше (мы реже выполняем инструкции типа mulsd, addsd, и больше mulpd, addpd). 

Для меня неожиданным оказалось то, что при реализации инвертирования матриц пришлось совершенно по другому организовать вычисления- способом, о котором мне не рассказывали в университете на лекциях по линейной алгебре и выч-мату, и который сильно отличается от “компилируемого”- но возник он только после того, как я уперся в ограничение на объем памяти в 256 байт (!) сидя на машине с 10-ю ядрами, 20-ю потоками и 32-мя гигами оперативы. 

Стоит ли останавливаться на достигнутом? Сейчас да. Самая низкая скорость обработки данных, полученная нами- это ~3ГБ/с на обращении матриц 4х4. 10 ядерный Xeon 2660v2 на этой задаче может прожевать 25GB/s, что дает суммарный прокачиваемый через шину памяти поток данных в 50GB/s, в то время, как четырех-канальный контроллер памяти DDR3-1866 позволяет примерно столько-же- ~59GB/s в теории, 45-50 в реале. То есть, даже на самых тяжелых расчетах мы упираемся не в ЦП, а в память. Переход на более быструю память DDR4 приведет и к смене процессора на более быстрый и с большим количеством ядер, и как показала практика- Ryzen 3900X  на DDR4-3200 и Xeon 2678v3 - оба все равно продолжают обрабатывать данные быстрее, чем их успевает подавать ОЗУ. 1 поток на Ryzen 3900x инвертирует ~5GB/s матриц 4х4, но процессор имеет лишь двухканальную память DDR4-3200, что позволяет ему ворочать те же самые 50GB/s, то есть, на таких расчетах его можно загрузить почти под завязку- но именно что почти- даже с DDR4-3600 в тестах скорость копирования данных в ОЗУ составляет ~56GB/s- то есть, как минимум одно ядро из 12 будет простаивать, вот и получается, что для нашей векторной алгебры, нужной для метода конечных элементов и аппроксимации частных производных методом наименьших квадратов- дальнейшие оптимизации или переход на AVX не приведут к повышению быстродействия. И даже переход на 4-х канальную DDR4-2400 не меняет ситуацию. А если нет разницы- зачем платить больше? Отсюда же понятно, что турбобуст мне без надобности.

Статья немножко опухла, и в нее не влезли векторные произведения векторов, быстрое решение СЛАУ 3х3 и 4х4, умножения разреженных матриц и еще под сотню велосипедных запчастей, но базовые примеры и самое страшное- инвертирование матриц- я описал.