SSE

Ну вот, наконец-то добрался до SSE. И предварительные результаты несколько обескураживающие. Я, конечно не ожидал, что использование векторных команд увеличит производительность пропорционально размеру вектора, все же даже решение СЛАУ, которое прямо просятся решаться векторными операциями, на 100% не векторизуется.

Но все же я надеялся хотя бы на двукратный рост на типе single (float), которому на SSE соответствует вектор из 4 элементов. На практике же максимальное ускорение относительно решения на Scalar SSE составило всего 1.66 раза, что на AMD FX-4350 соответствует решению 158 уравнений в секунду.
Очевидно, что на маленьких СЛАУ вышло даже небольшое замедление. На СЛАУ из 3 уравнений ускорение получилось 0.93.

Но самая печаль на типе double. В SSE ему соотвествует вектор всего из двух элементов. Максимальное ускорение достигнуто на СЛАУ из 192 переменных и составляет 1.18 раза. На маленьких размерностях (3 и 6 уравнений) ускорение меньше единицы, далее растет, после максимума на 192 уравнениях плавно снижается, приближаясь к единице 1536 уравнениях.

Одна из возможных причин состоит в том, что SSE не умеет выполнять арифметические операции с невыровненными на границу 16 байт данными. А при решении СЛАУ они гарантированно не выровнены.
AMD рекомендует в таких случая выполнять сохранение в две команды, записывая за раз по половине векторного регистра. Попробовал, ощутимой разницы не увидел: на каких-то размерностях стало быстрее, на других же наоборот.
Вроде бы в AVX эту "недоработку" пофиксили и арифметические команды выполняют и над невыровненными данными.
Тут кстати обнаружил, что моя версия Delphi не умеет в AVX. Надо подумать, как красивее обойти это органичение.

Ну и напоследок пример кода для double, возможно, мне опять какой-нибудь очевидный косяк замылился. 

procedure MakeDoubleNextRowsP_SSE;
asm
// R8 - P, R9w - n, R10w - index
  lea EBX, [R10d+1]; // BX - i
  cmp BX, R9w;
  jae @exit;
@Loop1:
  mov RCX, [R8 + RBX*8]; // RCX - S
  lea EDX, [R10d+1];
  lea EAX, [R9d-7];
  movsd XMM0, [RCX+R10*8];
  unpcklpd XMM0, XMM0;
  cmp EDX, EAX;
  jg @Loop2_4;
@Loop2:
  movapd XMM1, XMM0;
  movupd XMM2, [R11+RDX*8];
  mulpd XMM1, XMM2;
  movupd XMM2, [RCX+RDX*8];
  subpd XMM2, XMM1;
  movupd [RCX+RDX*8], XMM2;
  movapd XMM1, XMM0;
  movupd XMM2, [R11+RDX*8+16];
  mulpd XMM1, XMM2;
  movupd XMM2, [RCX+RDX*8+16];
  subpd XMM2, XMM1;
  movupd [RCX+RDX*8+16], XMM2;
  movapd XMM1, XMM0;
  movupd XMM2, [R11+RDX*8+32];
  mulpd XMM1, XMM2;
  movupd XMM2, [RCX+RDX*8+32];
  subpd XMM2, XMM1;
  movupd [RCX+RDX*8+32], XMM2;
  movapd XMM1, XMM0;
  movupd XMM2, [R11+RDX*8+48];
  mulpd XMM1, XMM2;
  movupd XMM2, [RCX+RDX*8+48];
  subpd XMM2, XMM1;
  movupd [RCX+RDX*8+48], XMM2;
  add DX, 8;
  cmp EDX, EAX;
  jle @Loop2;
@Loop2_4:
  add EAX, 4;
  cmp EDX, EAX;
  jg @Loop2_2;
  movapd XMM1, XMM0;
  movupd XMM2, [R11+RDX*8];
  mulpd XMM1, XMM2;
  movupd XMM2, [RCX+RDX*8];
  subpd XMM2, XMM1;
  movupd [RCX+RDX*8], XMM2;
  movapd XMM1, XMM0;
  movupd XMM2, [R11+RDX*8+16];
  mulpd XMM1, XMM2;
  movupd XMM2, [RCX+RDX*8+16];
  subpd XMM2, XMM1;
  movupd [RCX+RDX*8+16], XMM2;
  add DX, 4;
@Loop2_2:
  add EAX, 2;
  cmp EDX, EAX;
  jg @Loop2_1;
  movapd XMM1, XMM0;
  movupd XMM2, [R11+RDX*8];
  mulpd XMM1, XMM2;
  movupd XMM2, [RCX+RDX*8];
  subpd XMM2, XMM1;
  movupd [RCX+RDX*8], XMM2;
  add DX, 2;
@loop2_1:
  cmp DX, R9w;
  ja @Loop2_Fin;
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8];
  movsd XMM2, [RCX+RDX*8];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8], XMM2;
@Loop2_Fin:
  inc BX;
  cmp BX, R9w;
  jb @Loop1;
@exit:
end;