Стероидный велосипед: векторная алгебра, на ассемблере, в Delphi

  • Tutorial

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

Реклама
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее

Комментарии 33

    0

    В очередной раз показано как в Delphi, предназначенном для быстрой разработки GUI-приложений с доступом к БД, можно делать абсолютно разные вещи, в том числе быстрые математические расчёты.


    Кстати, ассемблерные команды SSE с расширениями 4.2, завезли в одну из недавних версий Delphi. До этого приходилось писать ассемблерные вставки в опкодах. Команд AVX в Delphi до сих пор нет. https://m.habr.com/en/post/441392/

      0
      SSE4.2 в ней давненько- все выше описанное компилируется даже в XE4
        0
        Да, вы правы, не так и недавно. Кстати, в Delphi 10.5 (2-я половина 2021) запланированы Math Performance Improvements. Не знаю, что это, но надеюсь то что я думаю.
      0
      Замечательно! Среагировал на статью, так как сам на днях (точнее лунной ночью) обращал матрицу 4х4. Результат использовался для некоторого урматфиза.
      Юмор
      Оборотень: «Не подскажете, какая сегодня фаза луны?»
      Я: «Полнолуние».
      Оборотень: «Спасибо».
      Я: «Не за что, обращайтесь».
          0
          Зачем в первом примере xor'ы для r13-15, они же гарантировано перепишутся полностью?
          Кстати, я пробовал еще EasyCode c UASM64 в паре с Delphi, тоже весьма удобно, особенно для инструкций, которые Delphi asm не поддерживает. В Delphi линкуется obj файлик из UASM64 без проблем.

          Например в EasyCode пишем -
          OurASMObj.asm
          ...
          TestProc1 Proc FastCall Frame i:Byte
          Xor Rax, Rax
          Mov Al, i
          Xor Al, 0xFF
          Ret
          TestProc1 EndP



          И затем в дельфи -
          ....
          {$L OurASMObj.obj}
          function TestProc1(a: Byte): Byte; external;
          
          implementation
          
          {$R *.dfm}
          
          procedure TForm1.Button1Click(Sender: TObject);
          var res: byte;
          begin
            res := TestProc1($AA);
            Edit1.Text := inttostr(res);
          end;
          end.

            0
            если поддерживать сразу две версии- 32 и 64, то проще принудительно обнулять все используемое в полном объеме, но Вы правы, в данном случае это избыточная операция.
            Использовать внешние редакторы ассемблера- на мой взгляд дело вкуса. но для tutorial'а- это было бы через чур. Да и для себя я в этом не вижу смысла- IDE, редактор и дебагер меня вполне устраивают.
            0

            Не думали для работы с матрицами 3х3 и 4х4 использовать видеокарту? Вроде как они адаптированы для работы с подобными матрицами...

              0

              В матрице 4x4 всего 16 элементов, столько же, сколько у автора машина умеет за одну инструкцию. А вот для видеокарты мне кажется размер матрицы маловат (и эффективность будет низкая)

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

                  Ах да, у вас даблы… жирненько. А чем обусловлено использование double? Нет ли возможности хотя бы на float/single заменить?

                    +1
                    на fp32 расчет деформаций идет с точностью порядка 10-6- а это по уровню- соответствует термическим деформациям от нагрева на 2-5 градусов. если хотим рассчитывать термодеформации- то приходится сидеть на fp64.
                0
                пробовали. но матрица у меня fp64, видеокарта такое не любит, это раз, а во вторых- я в любом случае на десяти-двенадцати ядрах полностью утилизирую ПСП и еще остается- поэтому карта мне не дает ускорения, она тоже упирается в ПСП. чтобы иметь прирост с карты- мне нужно закинуть в нее пару гигов данных, долго-долго их на ней крутить, и забрать с нее потом конечный результат. а если я постоянно гоняю данные туда-сюда- то профита нет, наоборот- медленнее получается.
                  0

                  Согласен… Лет 15 назад воял на делфи 7 трехмерную игру. Как раз в то время появились шейдеры. Работа с ними осуществлялась через язык HLSL, исполняемый в том числе в делфи. В нем как раз можно было работать с матрицами 3х3 и 4х4. Во что это выросло сейчас — сказать не могу. Делфи в основном использую для управления лабораторными приборами и обработки поступающих с них данных.

                0
                Во-первых, операции с матрицами реализованы уже, наверное, тысячи раз и странно, что не нашлось ни одной подходящей библиотеки.
                Я сначала подумал про Intel IPP, но пишут, что для операций с мелкими матрицами надо использовать MKL. Эти библиотеки сейчас бесплатные. Можно даже найти не очень старые заголовки для Дельфи (2019 года).

                Во-вторых, если всё-таки хотите писать самостоятельно, рекомендую обратить внимание на компиляторы Си. Они могут генерировать эффективный ассемблер (в т.ч. с векторизацией) даже из самого примитивного кода безо всяких зачатков оптимизации. Вот например ваши формулы для расчёта обратной матрицы 3*3 обычным методом, просто «в лоб» скопированные в Си-код: gcc.godbolt.org/z/8c4b1G
                Получилось с виду неплохо, около половины команд в ассемблере векторная (c окончанием на «pd»), по длине примерно так же, как у вас. Если хотите, можем сравнить реальную скорость.
                  0
                  очень часто при оптимизации важно даже не сами команды а их расстановка под какой либо конкретный процессор вплоть до того что между командами вставляются пустые команды чтобы гармонизировать их загрузку в процессор.
                    +1
                    библиотеки-то как раз нашлись, сырцов не нашлось.
                    Предлагать перейти на компиляторы Си- это, простите, моветон, во-первых- это очевидная мысль, а во-вторых- можно, но переписывать проект ради оптимизации десятка-другого мелких функций, когда результат- не гарантирован- плохая затея. Что до M4x4- так тут вообще ситуация патовая- не разложив элементы по регистрам лично я вообще не догадывался, что обращение можно делать не вылезая за пределы исходной матрицы- а без этого- никакая векторизация компилятором не поможет, потому что он отлично векторизует отвратительный алгоритм- с посредственным итоговым результатом.
                    Они могут генерировать эффективный ассемблер (в т.ч. с векторизацией) даже из самого примитивного кода безо всяких зачатков оптимизации

                    так в этом-то и соль- примитивный код можно эффективно оптимизировать и векторизовать, но код не всегда примитивный- не зря же Intel придумала свои интринсики software.intel.com/sites/landingpage/IntrinsicsGuide/#
                    А в случае с обращением матрицы- я в принципе не могу родить код (ни на С, ни на Fortran, ни на Pascal), который можно упаковать в одни регистры- хоть как там будет компилятор оптимизировать его- у меня даже мысли не возникало кидать элементы результата вместо получившихся нулей под диагональю исходной матрицы, чтобы не лезть в стек и минимизировать операции.
                    Если хотите, можем сравнить реальную скорость

                    хочу. меня интересует скорость обращения массива из 1млн матриц 4*4*FP64 в 1 поток на 1 ядре Ryzen3900 или любого близкого конкурента- используйте любой компилятор, любой алгоритм обращения, любой язык.

                    ассемблерный код с gcc.godbolt.org/z/8c4b1G я пока тестирую, и че-то я удивлен- он у меня выдает 490МБ/с, в то время как прямой дельфишный- 3300.
                    пробежав по коду- у CLANG-а- 63 инструкции, у меня- 64 всего, (это с префетчами + 6 инструкций проверка обусловленности). Но главное- у меня только один divsd, а по Вашей ссылке- два дива, а divsd- очень медленная. просто ужасно медленная- в ней 40 тактов можно потерять как нефиг делать. ну и чтение данных- я использую movupd, а CLANG везде поставил movsd, + невыровненность данных режет скорость чтения с памяти. Вот и получается, что на таком простом коде компилятор сделал простой ассемблер, который, как Вы сказали, «с виду не плохо».
                      +1
                      библиотеки-то как раз нашлись, сырцов не нашлось.
                      А почему нельзя без сырцов?
                      Вы неявно используете массу dll из Винды, и для них тоже нет сырцов.
                      Предлагать перейти на компиляторы Си- это, простите, моветон
                      Ну да, дельфистам сразу слышатся отголоски холиваров «Дельфи vs С++».
                      Нет, я предлагаю менять не Дельфи, а ассемблер на Си.
                      А в случае с обращением матрицы- я в принципе не могу родить код (ни на С, ни на Fortran, ни на Pascal), который можно упаковать в одни регистры
                      В регистры AVX (или скажем AVX-512) влезает гораздо больше.
                      хочу. меня интересует скорость обращения массива из 1млн матриц 4*4*FP64 в 1 поток на 1 ядре
                      Поскольку это ваша задача, то давайте вы напишете тестовую программу на Дельфи, а я к ней прикручу сишный вызов. Вам виднее, какие матрицы должны быть, какие примеры данных. Это и как приложение к статье будет полезно.
                      Но главное- у меня только один divsd, а по Вашей ссылке- два дива
                      С делением нехорошо получилось, да. Но это же легко правится вручную, дописать одну строку D = 1 / D и поменять деление на умножение. Или переключиться на систему команд AVX (-mavx) или AVX2/FMA (-march=haswell), будет одно деление.

                      В целом я не исключаю, что может быть медленнее, в конце концов, я на эту «оптимизацию» потратил 2-3 минуты. А вы на свою сколько дней?
                      Даже если придётся векторизовать вручную, на Си это будет компактнее, читабельнее и более гибко, можно легко переключаться между 32/64 битами и относительно легко — между наборами команд.
                        0
                        Вы неявно используете массу dll из Винды, и для них тоже нет сырцов

                        потому-что дальний вызов. а я где-то выше написал, что просто за вызов у меня штраф, да, я использую виндовые dll- но где? что-то редкое вывести, поток запустить-остановить, раз в пол-часа- данные на диск скинуть, то есть- очень редкие операции- потрачу я по две лишних секунды раз в пять минут- и не замечу. А векторы я свои делю на матрицы- миллиарды раз во все доступные потоки- и на этих миллиардах каждая мкс задержки- накапливается вполне ощутимо.
                        Еще потому-что сильно не хочется завязываться на чьи-то лицензии- сейчас нас стали проверять на чистоту всего используемого.
                        Ну и третье- очень было интересно на личном опыте пощупать: с одной стороны, общее мнение «в интернете»- что компиляторы стали такие умные, что оптимизируют лучше любого программиста, а с другой стороны- периодически выходят скромные статьи с разбором профайлингов и сказками про то, как какие-то одинокие самоучки делают умножение матриц на CUDE лучше, чем отдел разработки NVidia. На хабре было несколько примеров: умножение больших матриц, ракетный велосипед для преобразования FP64.toString.

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

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

                        заметано. за выходные выложу сюда линк на гитхаб для игрища.
                          0
                            0
                            github.com/RCherepanov/SSE_for_HABR наверное? Нет, пустой.
                            Можно и просто архивом выложить, без Гитхабов.
                              0
                              ой. я его по ошибке приватным сделал, теперь публичный
                              github.com/RCherepanov/SSE_bench_habr
                            0
                            периодически выходят скромные статьи с разбором профайлингов и сказками про то, как какие-то одинокие самоучки делают умножение матриц на CUDE лучше, чем отдел разработки NVidia. На хабре было несколько примеров
                            Мне кажется, конкретно операции с матрицами — очень ходовая и часто используемая штука, код стандартных библиотек должен быть вылизан до блестящего состояния и соревноваться с ними сложно.
                            Даже Ермолаев, при всей продвинутости его оптимизации, пишет в комментариях, что предположительно уступает MKL 5-10%.
                            Возможно, в вашем случае действительно будут влиять накладные расходы на вызов, у вас матрицы мелкие. В MKL для них предусмотрен какой-то хитрый инлайн, но при вызове из Дельфи он конечно работать не будет. Хотя можно написать промежуточную dll на Си, которая реализует функцию обработки массива матриц, с заинлайненной ф-ей MKL.
                            Мне также попадался одиночка, который заявляет, что у него быстрее, но для матриц общего вида — незначительно, те же 10% (transform inverse это, как я понимаю, упрощённый алгоритм для матриц трансформации в 3D-графике)
                            про читабельность на Си- я смотрел всякие портянки из интринсиков- честно- не вижу я там читабельности
                            Интринсики не особо читабельны, согласен.
                            Я обычно либо пытаюсь подтюнить код под автовекторизацию (это возможно, если понимать её логику и ограничения), либо использую векторные расширения Clang/GCC, они позволяют писать в «шейдерном» стиле. Ну знаете, в шейдерах/OpenCL/CUDA есть типы вроде float4, с которыми можно делать математику как с обычными float. Здесь тот же принцип.
                            Было
                                    __m128 sss = _mm_set1_ps(0.5);
                                    for (x = xmin; x < xmax; x++) {
                                        __m128i pix = _mm_cvtepu8_epi32(*(__m128i *) &lineIn[x]);
                                        __m128 mmk = _mm_set1_ps(k[x - xmin]);
                                        __m128 mul = _mm_mul_ps(_mm_cvtepi32_ps(pix), mmk);
                                        sss = _mm_add_ps(sss, mul);
                                    }

                            Стало
                                    float4 sss = (0.5);
                                    for (x = xmin; x < xmax; x++) {
                                        float4 pix = __builtin_convertvector(lineIn[x], float4);
                                        sss += pix * k[x - xmin];
                                    }

                            Кстати, подобную штуку пытаются сделать разработчики FPC, но у них она практически не работает, только на элементарных примерах вроде a=b+c.
                            0
                            Я погонял на выходных код, который Вы прислали (ассемблер вставил в дельфи).
                            Никак не могу заставить делфи класть передаваемый аргумент в RSP. Я не доконца понял, но судя по ассемблерному коду, CLang складывает копию матрицы в стек, а потом из стека забирает их снова- на этом должен быть штраф, но какой конкретно- не могу сейчас оценить. Дельфи упорно передает через регистр адрес первого элемента матрицы, и внутри выдергивает данные напрямую из памяти. Когда я Ваш пример заставил работать с RCX- он стал вполне себе быстрый:

                            «M3: Invert.compiler»: CalcTime =227263 mks; throutput 3543,5 MB/s
                            «M3: Invert_delphi_assembler_cleared»: CalcTime =196167 mks; throutput 4105,2 MB/s
                            «M3: Invert_gauss»: CalcTime =240986 mks; throutput 3341,7 MB/s
                            «M3: T_SSE.Invert_gauss»: CalcTime =182678 mks; throutput 4408,3 MB/s
                            «M3: T_SSE.Invert.direct»: CalcTime =135077 mks; throutput 5961,8 MB/s
                            «M3: T_SSE.Invert(N)»: CalcTime =145853 mks; throutput 5521,4 MB/s


                            «M3: Invert(packed matrix M3x3, compiler)»: CalcTime =195586 mks; throutput 3088,1 MB/s
                            «M3: T_SSE.Invert_CLANG»: CalcTime =161272 mks; throutput 3745,1 MB/s

                            Пока мои выводы такие:
                            a. (из разряда- удивительное рядом) компилятор делфи на простых математиках неплох.
                            б. если из ассемблера, сгенерированного делфи, выкинуть лишние операции (он там накидвает каких-то бесполезных пересылок)- то получается еще лучше.
                            в. ручная оптимизация с SSE- даже на простой математике дает заметный прирост, и получается лучше, чем у компилятора.
                            г. СLANG- дает вполне компактный и эффективный код, но все равно может пропустить очевидные вещи (как с повторным divsd на одних и тех же данных). между CLang v5.0 и CLang v11.0.1- разница видна невооруженным взглядом.
                            д. CLang v11.0.1 лучше ICC- то, чего нагенерил ICC- вообще никому показывать нельзя- там 5 divpd/sd вместо 1го!
                            е. выровненность данных дает небольшой прирост. M3: «M3: Invert.compiler» vs Invert(packed matrix M3x3, compiler)"

                              0
                              г. СLANG- дает вполне компактный и эффективный код, но все равно может пропустить очевидные вещи (как с повторным divsd на одних и тех же данных)
                              div — это так, досадное недоразумение. Не в div-ах суть. Я даже не буду спрашивать, сколько div-ов на таком коде генерирует Дельфи :) Нормальный разработчик делает замену в коде множественного деления на умножение «на автопилоте».
                              А суть вот в чём: я переделал код на обработку массива матриц, и теперь векторизатор отработал лучше. Он считает по 2 матрицы за итерацию, за счёт этого все команды цикла векторные.
                              Причём можно сделать загрузку матриц в регистры эффективнее, если хранить их по-другому: не как массив структур, а как раздельные массивы для каждого элемента матрицы. Смотрите сами, насколько это допустимо в вашей ситуации.
                              Раздельные массивы хороши тем, что лучше масштабируются на любую систему команд, можно считать по 4 матрицы за итерацию с AVX.
                              И всё это в 20-30 строчках максимально простого кода, никаких ассемблерных портянок вручную.
                                0
                                #pragma clang loop vectorize(assume_safety)
                                мне такое очень не желательно, я часто сохраняю инвертированную матрицу поверх исходной.
                                Я согласен с Вами в том смысле, что можно сделать простой код, который будет хорошо компилироваться специальным компилятором (CLang 11.0.1, но не GCC10.2 например, и не msvc или icc, я до общения с Вами пребывал в уверенности, что лучшая оптимизация- в ICC- как никак, компилятор от производителя процессора! а оно вот оно че оказывается). Вполне себе оправданный подход в определенных условиях (например, когда можно ради упрощения жизни компилеру переломать базовые структуры данных- ведь на эту самую Mat3x3 завязана половина кода, ее тронь- там столько всего полезет!). Но совершенно спокойно можно прямо в родной старючей delphi xe4 (не слезая с кактуса) на старючем SSE4.2 без AVX забить всю ПСП и пережевывать данные быстрее, чем они через нее пролазят. КМК, туториалы для того и нужны, чтобы показывать разные способы решения одних и тех же задач, чтоб можно было выбирать под свои нужды. Лично для меня одной из причин залезть именно в ассемблер было то, что предыдущий опыт использования «ускоряющей» dll был полон мучительной отладки, когда какие-то данные вызвали ошибку в потрохах ДЛЛ, а дебагер не может до них докопаться. А когда залез- то оказалось, что это совсем не страшно. :-)
                                  0
                                  до общения с Вами пребывал в уверенности, что лучшая оптимизация- в ICC- как никак, компилятор от производителя процессора!
                                  ICC тоже хорош, но у него своя система параметров, и не все из них работают через godbolt. Попробуйте -fast, и будет таки fast, но принудительно включится FMA/AVX2. В общем, я толком не знаю, как с ним обращаться.
                                  когда можно ради упрощения жизни компилеру переломать базовые структуры данных
                                  Да всё уже, отбой, не надо ломать. На практике раздельные массивы «не взлетели». Наверное, виноват не последовательный доступ к памяти, лезем в память по 9-и разным указателям вместо 1-го. Может быть, если бы хранить элементы мини-массивами по 4 штуки, было бы лучше, но это совсем уже неудобно в работе.
                                  От новых инструкций (AVX, FMA) тоже толку немного.
                                  Но в результате я всё-таки обогнал ваш ассемблер примерно на 30%, а для компактной матрицы сильно, раза в 3.
                                  Аккаунта на Гитхабе нет, поэтому архивом, там же и результаты.
                                  забить всю ПСП и пережевывать данные быстрее, чем они через нее пролазят
                                  Не уверен, что если вы забиваете ПСП, то это повод для гордости. Может слишком мало арифметики на чтение/запись? Но да-а, теперь уже не хочется ничего менять, когда всё захардкожено ассемблером.
                                  Лично для меня одной из причин залезть именно в ассемблер было то, что предыдущий опыт использования «ускоряющей» dll был полон мучительной отладки, когда какие-то данные вызвали ошибку в потрохах ДЛЛ, а дебагер не может до них докопаться.
                                  Если это ваша dll-ка, то отлаживать её можно из сишной IDE c дельфийским host application.
                                  Это добавляет сложности, но и отладка ассемблера — тоже так себе удовольствие, особенно через пару лет после того, как вы этот ассемблер писали.
                                  Но совершенно спокойно можно прямо в родной старючей delphi xe4 (не слезая с кактуса)
                                  Хорошо, не слезайте, больше не буду уговаривать :)
                                    0
                                    коллега! не уходите так рано! помните- у меня главная засада была не в 3х3, а в 4х4!
                                      0
                                      У вас для Гаусса 4x4 нет готового кода на Дельфи, который можно взять за образец, либо ассемблер, либо словесное описание.
                                      Здесь я взял за основу 3x3, может сами допишете?
                                      gcc.godbolt.org/z/1KKG9a
                        0
                        Скажите пожалуйста сколько времени потратили на переписывание итого?
                        Покрывали ли как либо тестами? Как я понимаю по хорошему надо тестами и проверять что корректно восстанавливаются регистры при выходе и не портят не используемые данные?
                        Идентичны ли результаты полученным на обычном фпу?

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

                        Разве щас не под капотом дельфи компилятор llvm?
                          0
                          чистого времени на все SSE было потрачено около двух недель- это порядка 70-ти различных функций. Но большая часть из них очень похожа (типа, a.i += b.ij*cj and a.i += b.ij*c.j*k)- они копипастились с минимальными правками, и основная возня была именно с матрицей 4*4. Корректность расчетов покрыта unit-тестами, идентичность результатов с fpu- 50/50, где-то побитовая точность, где-то- в пределах погрешности округления, ибо перемена порядка действий влияет на последние знаки результата.
                          Про подкапотное пространство не скажу- на работе лицензия XE4, там все грустно, llvm еще даже в проектах не было, и современный FPC дает существенно более быстрый код. Что щас в Берлине- не смотрел, на ноуте он стоит, но чисто символически- как запасной аэродром, тестить на нем производительность даже не пробовал (N4200 не для расчетов никак).
                            +1
                            Да сейчас уже не Берлин давно… Уже несколько лет не Берлин (10.1). Сейчас Сидней (10.4).
                            +2
                            LLVM для ARM (iOS, Android) и Linux. Для Windows — старый много лет допиливаемый компилятор.

                          Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                          Самое читаемое