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