29.09.2020

Продолжаю оптимизировать

Точнее, пытаюсь закончить оптимизировать Scalar SSE в тех местах, которые до конца не проверил, или проверил не правильно. 

1. Замена вычитания сложением. В чем суть этой попытки оптимизации? У нас есть вот такое выражение:
  S[j] := S[j]-R[j]*S[index];
На ассемблере я его транслировал в такой набор команд
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8];
  movsd XMM2, [RCX+RDX*8];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8], XMM2;

в двух операторной форме команд, принятой в x86, я вынужден сначала загрузить S[j] в регистр, а уже потом выполнять вычитание.
Если же у S[index] перед внутренним циклом поменять знак, то операцию вычитания можно заменить сложением без предварительной загрузки из памяти в регистр XMM2. В результате получится вот такой код:
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8];
  addsd XMM1, [RCX+RDX*8];
  movsd [RCX+RDX*8], XMM1;
Как видим, он на одну команду короче, и, возможно, быстрее. Не стоит забывать, что нам пришлось во внешней цикле добавить команду изменения знака S[index].
Но будет ли этот код реально быстрее? Зависит от того, какая команда выполняется быстрее: умножение XMM1 на значение в память или загрузка в XMM2 значения из памяти. Дело в том, что команда addsd не может начаться, пока не закончится умножение XMM1.
Поэтому, если умножение медленное, то XMM2 успеет загрузиться из памяти (так как эти команды независимы) и затем можно будет вычесть из XMM2 XMM1 без потери производительности. То есть в этом случае замена вычитания сложением не нужна, так как ничего не даст.
В том случае, если умножение выполняется быстрее, чем загрузка XMM2, то ситуация будет обратная: выгоднее будет второй вариант со сложением.
На мой взгляд, разницы тут не должно быть много, но давайте посмотрим на результаты.
На типе single (float) ускорение от замены вычитания на сложение составило от сотых долей процента на системах из 3 уравнений, до 4-5% на больших системах до 48 переменных включительно. А вот на больших размерностях результат оказался парадоксальный: незначительное снижение производительности на системах из 96 уравнений и потери до 3.5% и выше на больших размерностях.
На типе double ускорение не превышает 3%, достигая максимума на 24 уравнениях, затем постепенно снижается и заменяется замедлением на системах в 1536 уравнений.
В общем, я не увидел особого смысла в такой оптимизации.

2. Развертывание цикла. Я прошлый раз ограничился четырьмя операциями, но Иван предложил попробовать увеличить их до восьми.
Результат моих размышлений на эту тему.
Ускорение от развертывания можно получить за счет снижения доли служебных операций на организацию цикла. Но это работает только в том случае, когда тело цикла (рабочая нагрузка) очень проста и по времени сравнима с работой по организации самого цикла.
В моем случае сам цикл организуется всего тремя командами: двумя короткими регистровыми операциями и условным переходом. А полезной нагрузкой является 5 команд, в том числе три с памятью. В принципе, доля служебных команд достаточна велика и такой цикл имеет смысл разложить на две итерации. Дальнейшее разложение приведет лишь к незначительному росту производительности, тем более, что хорошо работает блок предсказаний переходов, команды организации цикла и полезной нагрузки исполняются на разных блоках, так что частично могут выполняться параллельно.
Второй источник ускорения связан с отсутствием сбросов конвейера. Чем длиннее линейная часть программы, без условных переходов, тем полнее заполнен конвейер и выше производительность процессора. Наличие условных переходов сбрасывает конвейер,что приводит к сильному снижению производительности.
Вышеописанная страшилка более актуальная для очень старых процессоров. Современные версии процессоров имеют очень хороший блок предсказания переходов, который работает совместно с процессором и остановка конвейера происходит только в том случае, если блок предсказаний ошибся. В частности, для for-циклов блок предсказаний ошибается очень редко, в идеальном случае лишь один раз ‒ когда цикл заканчивается. Думаю, что именно этот источник ускорения при решении моей задачи имеет очень низкое значение.
Третий источник ускорения
‒ суперскалярная архитектура, причем именно она, на мой взгляд, дает наибольший эффект в моем случае. Посмотрим на полезную нагрузку цикла. В случае одной операции на итерацию у нас требуется выполнить одно умножение и одно сложение, причем сложение зависит от умножение. Фактически здесь нет операций, которые могли бы выполнятся параллельно.
Если выполнять две операции за итерацию, то нужно выполнить 2 умножения и 2 сложения, тут уже есть что распараллелить. 2 умножения можно выполнять параллельно, но со сложениям из-за зависимости от умножений ситуация сложнее, скорее всего параллельно умножениям удастся выполнить лишь одну команду сложения.
При увеличении количества операций на итерацию ситуация с распараллеливанием становится все лучше и лучше. Теоретически могут выполняться 4 инструкции за такт, но практически, есть подозрение, их количество лежит в пределах от 2 до 3, причем ближе к двум, по всей видимости. Это связано с тем, что умножение все же выполняется медленнее сложения, поэтому для сумматоров данные всегда поставляются с задержкой.
Тем не менее, попробуем оценить, насколько увеличивается производительность при переходе от 4 операций на итерацию к 8.
На типе single (float) ускорение составило от 0% (3 переменных) до 11% (192 переменных), затем уменьшилось до 4% на 384 переменных.
На типе double ускорение максимум составило около 7% на 192 переменных, на всех остальных размерах ускорение составляет около 1%, что ниже погрешности измерений.
В результате решил оставить развертывание на восемь операций для типа
single (float), а для типа double остановиться на четырех.

3. Решил побаловаться с командами PREFETCH. Эффект для моей задачи равен нулю, так как, видимо, механизм автоматической предзагрузки кэша в этом случае справляется на ура.

27.09.2020

Оптимизирую Scalar SSE

Первый вариант моей оптимизации путем развертывания цикла оказался не очень удачным, как совершенно справедливо указал Иван Колесников в комментариях к посту Scalar SSE. Так что исправляю.

1. Да, организация цикла были весьма корява, исправление существенно увеличило производительность.

2. Ручное переименование регистров, по всей видимости, сбивало с толку механизм, встроенный в процессор. Отказ от него также привел к увеличению производительности.

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

4. А вот предложение заменить вычитание умножением на -1 и сложением наоборот, замедляло производительность.

5. Не стал реализовать концовку аналогом оператора case (switch). Вот мои соображения:
во-первых, такой оператор не очень быстр: тут либо несколько условных переходов, либо, более быстро косвенный безусловный переход по адресу в памяти (8-ми байтному); не уверен, что адрес будет в кэше, предыдущий длинный цикл по строке скорее всего его вытеснит; и переход такой будет лишь один раз на все итерации внутреннего цикла;
во-вторых, такой вариант не очень хорошо ложится на последующий переход к векторному варианту решения.
Спорно, конечно, но решил что этот вариант не сильно ускорит решение СЛАУ.

Насколько удалось ускорить, покажут будущие тесты.

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

Ну и напоследок код главного цикла (цикл для верхней строки тоже аналогично оптимизировал; ну и то же самое сделал для FPU):

procedure MakeDoubleNextRowsPASM;
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];
  mov EAX, R9d;  // --
  movsd XMM0, [RCX+R10*8];
  sub EAX, 3;
  cmp EDX, EAX; // --
  jg @Loop2_2;
@Loop2:
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8];
  movsd XMM2, [RCX+RDX*8];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8], XMM2;
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8+8];
  movsd XMM2, [RCX+RDX*8+8];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8+8], XMM2;
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8+16];
  movsd XMM2, [RCX+RDX*8+16];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8+16], XMM2;
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8+24];
  movsd XMM2, [RCX+RDX*8+24];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8+24], XMM2;
  add DX, 4;
  cmp EDX, EAX; // --
  jle @Loop2;
@Loop2_2:
  add EAX, 2;
  cmp EDX, EAX;
  jg @Loop2_1;
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8];
  movsd XMM2, [RCX+RDX*8];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8], XMM2;
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8+8];
  movsd XMM2, [RCX+RDX*8+8];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8+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;

24.09.2020

Scalar SSE. Предварительные результаты.

Оттестировал разные процессоры на пока еще не до конца оптимизированном коде. Думаю, после оптимизации он станет существенно быстрее. Но так как оптимизация будет инвариантна для всех процессоров, то относительный результат не сильно изменится.

Во-первых, развертывание циклов даже не в самом оптимальном варианте показало все же более серьезное ускорение, особенно на больших размерностях.
У AMD FX-4350 это порядка 30% для single (float) и 36% для double. Мобильный AMD A10-4600M показывает чуть худший результат: 20% и 33% соответственно. Рекордсмен среди процессоров от AMD, конечно Ryzen 5, у него 41% и 67%. Только надо учитывать, что это относительно производительности не развернутого цикла FPU (x87), который в Ryzen, как выяснил Иван (см. комментарии к прошлому посту), зарезан в 2 раза по отношению к Scalar SSE, тут все четко получилось.
Intel i3-3227U показывает 29% и 27%. Старшие Intel тут менее результативны: Intel i7-6700HQ 20% для обоих типов и Intel i5-8300H 18% и 17%.

Теперь сравнение производительности Scalar SSE для разных процессоров относительно AMD FX-4350.

 

Как видим, появился тест, в котором AMD Ryzen 5 3600 смог догнать интеловские процессоры. Правда, радость омрачает тот факт, что догнать он смог только мобильные процессоры.
Люблю я процессоры AMD, в первую очередь, конечно, за цены. Еще они первыми добавили поддержку PCIe 4.0, что тоже плюс, правда, не очень пока ясно с чем его кушать, этот плюс.
Но справедливости ради должен сказать, что по однопоточной производительности они все еще далеки до отчетливого лидерства. Но надо проверить и другие дисциплины.😉

Еще бросается в глаза тот факт, что на младших мобильных, что Intel, что AMD, в многозадачном режиме SSE приходится жить похуже, а вот старшие мобильные Intel гораздо лучше с этим справляются. Даже лучше, чем с FPU. На Ryzen многопоточный режим одинаково эффективен что на FPU, что на Scalar SSE, ну и банально: этот процессор однозначный лидер в многопоточной производительности тупо из-за числа ядер.

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

20.09.2020

Scalar SSE

Начал реализовывать. В принципе, на что рассчитывал? FPU (x87) осуществляет операции, расширяя предварительно операнды в памяти до 80 бит, в то время как Scalar SSE этого вроде бы не делает, выполняя операции точно в размер типа.
Соответственно, на типе single (float) в 32 бита теоретически можем получить ускорение в 2.5 раза. Заменяем команды FPU на соответствующие команды Scalar SSE и получаем... ту же самую производительность, по крайней мере на
AMD FX-4350. То есть производительность FPU и Scalar SSE отличается на плюс/минус погрешность измерения.
Хм... Крайне странно. Попробую что ли развернуть циклы?

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

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

В принципе, возможно и другое аппаратное решение: одно АЛУ, которое выполняет действие над всем регистром SSE, в результате которого получается сразу четыре результата (для типа single). В этом случае раскрытие цикла на Scalar SSE не даст слишком большого выигрыша.

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

В общем, такое ощущение, что Scalar SSE ничего не знает про суперскалярность! 😉 Или я просто не умею его готовить?
Ниже пример кода, реализующего внутренний цикл, на этот раз, правда, для типа double. В принципе, ничем от single не отличается, кроме размера операндов. Результаты тоже аналогичны.

procedure MakeDoubleNextRowsPASM;
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];
  movsd XMM0, [RCX+R10*8];
  mov ax, R9w;  // --
  sub ax, R10w; // --
@Loop2:
  cmp ax, 4;
  jb @Loop2_2;
  movsd XMM1, XMM0;
  movsd XMM2, [RCX+RDX*8];
  movsd XMM3, XMM0;
  mulsd XMM1, [R11+RDX*8];
  movsd XMM5, XMM0;
  movsd XMM4, [RCX+RDX*8+8];
  subsd XMM2, XMM1;
  mulsd XMM3, [R11+RDX*8+8];
  movsd XMM7, XMM0;
  movsd XMM6, [RCX+RDX*8+16];
  movsd [RCX+RDX*8], XMM2;
  subsd XMM4, XMM3;
  mulsd XMM5, [R11+RDX*8+16];
  movsd [RCX+RDX*8+8], XMM4;
  subsd XMM6, XMM5;
  mulsd XMM7, [R11+RDX*8+24];
  movsd XMM8, [RCX+RDX*8+24];
  movsd [RCX+RDX*8+16], XMM6;
  subsd XMM8, XMM7;
  sub ax, 4
  movsd [RCX+RDX*8+24], XMM8;
  add dx, 4;
  jmp @Loop2_Fin;
@Loop2_2:
  cmp ax, 2;
  jb @Loop2_1;
  movsd XMM1, XMM0;
  movsd XMM2, [RCX+RDX*8];
  movsd XMM3, XMM0;
  mulsd XMM1, [R11+RDX*8];
  movsd XMM4, [RCX+RDX*8+8];
  subsd XMM2, XMM1;
  mulsd XMM3, [R11+RDX*8+8];
  movsd [RCX+RDX*8], XMM2;
  subsd XMM4, XMM3;
  sub ax, 2;
  movsd [RCX+RDX*8+8], XMM4;
  add dx, 2;
  jmp @Loop2_Fin;
@loop2_1:
  movapd XMM1, XMM0;
  mulsd XMM1, [R11+RDX*8];
  movsd XMM2, [RCX+RDX*8];
  subsd XMM2, XMM1;
  movsd [RCX+RDX*8], XMM2;
  inc DX; //--
@Loop2_Fin:
  cmp DX, R9W;
  jbe @Loop2;
  inc BX;
  cmp BX, R9w;
  jb @Loop1;
@exit:
end;


17.09.2020

Многопоточное ускорение

Небольшой анализ, насколько хорошо ускоряется решение множества СЛАУ в параллельном режиме. Сначала график.


Мелкая деталь, которая обращает на себя внимание: AMD FX-4350 и A10-4600M обеспечивают почти одинаковое ускорение с небольшим преимуществом десктопного процессора. Что не удивительно ‒ в основе лежит одинаковая архитектура PileDriver. Видимо, небольшой проигрыш определяется отсутствием кэша третьего уровня у мобильной версии.
Тем не менее, AMD здесь постаралась на славу, если вспомнить, что на 4 ядра в этой архитектуре имеется всего два полноценных АЛУ с плавающей точкой. Несмотря на это, в параллельном режиме обеспечивается трехкратное ускорение!

Но еще более впечатляющих результатов AMD добилась в AMD Ryzen 5 3600. Несмотря на 6 физических ядер, на небольших размерностях ускорение превышает 6 раз. Старшие мобильные Intel вроде тоже так умеют, но на более меньших размерностях, которые в этот раз в тест не попали.

Есть и общая закономерность, которую демонстрируют все процессоры: замедление ускорения от параллельного решения СЛАУ при росте размерности задачи. Особенно это заметно у десктопного Ryzen 5 3600, а также мобильных Intel i5-8300H и Intel i7-6700HQ. Причем последний начинает сдуваться даже раньше, на 384 переменных.
В итоге ускорение параллельной обработки у мобильных процессоров Intel становится меньше единицы, т. е. один поток выполнит работу быстрее, чем восемь!
Десктопные процессоры не демонстрируют такую особенность, видимо, это произойдет на одной из следующих размерностях, которые я не тестировал.

Объясняется эта закономерность несоответствием скорости процессора и пропускной способностью памяти. При решении СЛАУ каждая операция слишком проста, и когда строки памяти становятся слишком длинные, то каждое ядро слишком часто обращается к памяти, увеличивается промахи мимо кэша, а пропускной способности шины памяти становится недостаточно, что бы успеть обслужить запросы всех ядер.
Видимо, при решение подобных задач, обрабатывающих простым алгоритмом огромный объем памяти, использование технологии HyperThreading выглядит не очень оправданным. Т. е. ее видимо, либо необходимо отключить, либо просто использовать в два раза меньше вычислительных потоков, чем ядер.
Впрочем, это не гарантирует, что на очень больших размерностях этого будет достаточно. Все зависит от скорости шины памяти. Если она достаточна медленна, а обработка одного элемента слишком быстра, то количество потоков придется еще больше ограничить.

Кстати, самые слабенькие процессоры из теста, Intel i3-3227U и AMD A10-4600M демонстрируют крайне слабое падение ускорения в зависимости от размерности задачи. Они просто достаточно медленные, и ядер не так много, поэтому подсистема памяти успевает обслуживать их запросы.

15.09.2020

Новые результаты

Оттестировал новой 64-битной версией программки FPU (x87) на нескольких компьютерах.
Попутно понял, что не совсем правильно построил логику работы теста. Сейчас у меня, если памяти не хватает, то тест на этой размерности пропускается и берется, следующая, большая.
Логика тут простая: при увеличении размерности в 2 раза требуемый объем памяти увеличивается в 4 раза, а производительность падает (теоретически) в 8 раз, то есть при увеличении размерности нужно меньше СЛАУ, что бы провести тест.
Если же использовать только доступную память, то время решения будет слишком маленьким, а погрешность теста ‒ высокая.

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

Параллельно этим размышлениям появилась следующая мысль: видимо, при покупке нового компьютера, если есть цель обеспечить максимальную производительность современного процессора, нужно много памяти. Иначе процессор будет простаивать, ожидая загрузки страниц виртуальной памяти с внешнего ЗУ. Конечно, я тут не имею в виду процессорные обрубки для офисных пиш. машинок.
Минимум сейчас требуется, кмк, 16 ГиБ, рекомендуется 32 ГиБ, ну а если брать топовые процессоры с количеством ядер 8 и больше, то тут нужно все 64 ГиБ.
Для рабочих станции все цифры, мне кажется, надо увеличить минимум вдвое.

А теперь итоги теста.

Как видно, старенький мобильный i3-3227U с 2 ядрами выступает вровень с AMDшным A10-4600M с четырьмя ядрами. Более новый мобильный i5-8300H обгоняет i7-6700HQ какого-то более старого поколения.
А вот свежий бюджетный AMD Ryzen 5 3600 откровенно разочаровывает, так как новый десктопный процессор проигрывает обоим более пожилым мобильным конкурентам от Intel. Всегда почему-то думал, что десктопный процессор должен быть быстрее аналогичного по позиционированию мобильного.

После массы новостей, заполонивших тематические сайты, у меня создалось впечатление, что все же Ryzen должен выступать как минимум вровень с процессорами от Intel. Но нет. По крайней мере по FPU он явно проигрывает.
Правда, есть робкая надежда, что AMDшные процессоры хороши в каких-то иных дисциплинах. Проверю это чуть позже. Ну и понятно, что они точно дешевле Intel, а кроме того, у них есть PCIe 4.0.

Теперь посмотрим на многоядерную производительность в параллельном режиме.


Вот здесь процессор AMD наконец-то продемонстрировал на что он способен. Конечно, он одержал уверенную победу за счет большего количества ядер (6 физических против 4 у мобильных i5 и i7). Надо сказать, что AMD в Ryzen существенно улучшила многозадачность. Более подробно я напишу об этом в отдельном посте.
Интересно посмотреть на аутсайдеров. У i3 2 физических ядра и 4 виртуальных, которые дает Hyper Threading. У A10 вроде как 4 физических ядра и никакого HT. Но сделан A10 по технологии Piledriver, в которой лишь одно FPU на каждые два ядра (кстати, как и у FX-4350). Поэтому эти два процессора, i3 и A10, фактически и выступают на равных в многопоточном тесте, AMD даже оказался чуть-чуть быстрее. А extended (long double), в котором у Intel традиционно преимущество (видимо, память быстрее работе при доступе к невыровненным данным) я в этот раз не тестирую.
 
Ну и наконец диаграмма интегральной производительности.

10.09.2020

64 бита

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

В процессе обнаружил пару любопытных моментов.
Во-первых, так как не придумал способа, как сгенерировать гарантированно совместную СЛАУ, то просто включил в код проверку на то, что очередной диагональный элемент существенно больше машинного нуля для соответствующего типа данных.
В принципе, так даже правильнее, так как в реальности сходимость решения надо проверять всегда. На производительности это практически никак не сказалось.
А потом сразу проверил, насколько точно решаются СЛАУ. И выяснил, что тип single (или float в С-подобных языках) дает достаточно большую ошибку для СЛАУ больших размерностей, в районе несколько сотен переменных, ближе к тысяче.
Критерием для меня служила ошибка, превышающая 5% относительно точного решения. И больше половины решенных СЛАУ больших размерностей имеют максимальную ошибку больше 5%. Таким образом, я пришел к выводу, что СЛАУ размерностью более 500 переменных нет смысла тестировать на типе single.
С типом double же подобных проблем не было, по крайне мере для тех размерностей, что я смог проверить: максимальная ошибка была крайне мала.

Во-вторых, начал реализовывать тип extended (long double) и получил на ровном месте некоторую странную ошибку. Через какое-то время вспомнил, что AMD рекомендовало избегать x87 инструкций в 64-битном коде, заменяя их скалярными аналогами из SSE. А Embarcodero последовал этому совету весьма простым способом, элементарно приравняв тип extended к типу double в 64-разрядном режиме.
То есть код-то для решения таким способом у меня есть, а вот кода, который генерирует СЛАУ, получается нет, так как он на Pascal у меня написан. Сначала хотел описать свой тип, что бы таки реализовать тест на расширенном вещественным, но потом что-то заленился. Может быть, потом как-нибудь.

А зачем он нужен, этот расширенный тип? Банальный ответ: для повышения точности промежуточных вычислений. Но таких случаев в действительности не очень много. На практике, например, я разницу точности в решениях СЛАУ с использованием x87 и scalar SSE не увидел, хотя сильно глубоко там и не копал.

Что получилось в результате? Первое: 64-битный код оказался почти всегда быстрее 32-битного, причем на малых размерностях значительно (почти на 50%). В основном, оптимизация на большем количестве регистров позволила практически полностью отказаться от операций со временной памятью. Ну и хранение коэффициентов всех СЛАУ в одном массиве улучшило коэффициент попадания в кэш.
Но на больших размерностях все же 32-битный режим совсем чуть-чуть быстрее. Может, тут все же сыграли отрицательную роль 64-битные указатели?

Ну и выложу напоследок код, вдруг кто-то увидит недостатки и предложит улучшения.
Сначала типы, приведу только для систем на single:

  TSingleArray = array[0..65535] of single;
  PSingleArray = ^TSingleArray;
  PPointerArray = ^TPointerArray;
  TLESInfoP = packed record
    Data : Pointer;
    Systems : PPointerArray;
    Rows : PPointerArray;
    LESDataSize : UInt64;
    LESCount : UInt64;
    LESSize : word;
    ElementSize : byte;
  end;

TLESInfoP ‒ это структура, которая собирает в кучу общие сведения о том, сколько и какие системы хранятся в общем массиве Data.
Systems
‒ это массив указателей на Rows, i-ый элемент которого указывает на первую строку i-ой системы. Ну а Rows просто указывает строки.

На а теперь код, сначала на Pascal, который я затем перевел на ассемблер. Что бы понятнее было.

procedure SolveTestLES(var P : TLESInfoP; LESNum : UInt64);overload;
begin
  case P.ElementSize of
    4  : SolveSingleLESP(P.Systems[LESNum], P.LESSize);
    8  : SolveSingleLESP64(P.Systems[LEsNum], P.LESSize);
    10 : SolveSingleLESP80(P.Systems[LEsNum], P.LESSize);
  else
    raise Exception.Create('Wrong LES element size!');
  end;
end;

Это что бы было понятно, при чем тут TLESInfoP. На самом деле вся работа в SolveSingleLESP. Тут я допустил некоторую путаницу: Single в название, это типа одиночное уравнение решаем, а не используемый тип данных.

procedure SolveSingleLESP(P : PPointerArray; n : word);overload;
var
  i : word;
begin
  for i := 0 to n-1 do
  begin
    if FindSingleMaxP(P, i, n) <> 0 then
      exit;
    Make1RowP(P[i], i, n);
    MakeNextRowsP(P, i, n);
  end;
  for i := n-2 downto 0 do
    ReverseGJ(P, i, n);
end;

function FindSingleMaxP(P : PPointerArray; index, n : word):byte;
var
  i : word;
  Max : word;
  A : Pointer;
  x : extended;
begin
  Max := index;
  x := abs(PSingleArray(P[index])[index]);
  for i := index+1 to n-1 do
    if abs(PSingleArray(P[i])[index]) > x then
    begin
      Max := i;
      x := abs(PSingleArray(P[i])[index]);
    end;
  if Max <> index then
  begin
    A := P[index];
    P[index] := P[Max];
    P[Max] := A;
  end;
  if x < c_LowDiagonalSingle then
    Result := 1
  else
    Result := 0;
end;

procedure Make1RowP(P : PSingleArray; index, n : word);overload;
var
  i : word;
  d : extended;
begin
  d := 1/P[index];
  for i := index+1 to n do
    P[i] := P[i]*d;
end;

procedure MakeNextRowsP(P : PPointerArray; index, n : word);overload;
var
  i, j : word;
  S, R : PSingleArray;
begin
  R := P[index];
  for i := index+1 to n-1 do
  begin
    S := P[i];
    for j := index+1 to n do
      S[j] := S[j]-R[j]*S[index];
  end;
end;

procedure ReverseGJ(P : PPointerArray; index, n : word);overload;
var
  i : word;
  S : PSingleArray;
  d : extended;
begin
  d := PSingleArray(P[index+1])[n];
  for i := index downto 0 do
  begin
    S := P[i];
    S[n] := S[n] - d* S[index+1];
  end;
end;

Ну а теперь все тоже самое, только на ассемблере (логика немного изменена по сравнению с Pascal ‒ вместо общей процедуры, которая вызывает три конкретных, здесь просто три отдельных процедуры):

procedure SolveTestSingleLESASMP(var P : TLESInfoP; LESNum : UInt64); // P - RCX, LESNum - RDX;
asm
// Free modifiing RAX, RDX, RCX, R8, R9, R10, R11
// R8 - P, R9w - n, R10w - index
  push RBX;
  mov R8, P.Systems;
  movzx R9, word ptr P.LESSize; // в R9w - n;
  mov R8, [R8+LESNum*8]; // R8 указатель на указатель на первую строку СЛАУ RDX  теперь свободен
  xor R10, R10; // R10w - i;
  xor RBX, RBX; // готовим RBX для работы с индексами
@Loop1:
  call FindSingleMaxPASM;
  test al, al;
  jnz @exit;
  call Make1RowPASM;
  call MakeSingleNextRowsPASM
  inc R10w;
  cmp R10w, R9w;
  jb @Loop1;
  dec R10w;
@Loop2:
  dec R10w;
  call ReverseSingleGJPASM;
  test R10w, R10w;
  jnz @Loop2;
@exit:
  pop RBX;
end;

procedure FindSingleMaxPASM;
asm
// R8 - P, R9w - n, R10w - index
  fld dword ptr c_LowDiagonalSingle;
  lea RCX, [R10d+1];
  mov BX, R10w; // BX - Max
  mov R11, [R8 + R10*8];
  fld dword ptr [R11 + R10*4];
  fabs;
  cmp CX, R9w;
  jae @exit; // нечего вычислять
@Loop:
  mov R11, [R8+RCX*8];
  fld dword ptr [R11 + R10*4];
  fabs;
  fucomi ST(0), ST(1);;
  jbe @GL;
  fxch;
  mov bx, cx; // Max := i
@GL:
  fstp st(0);
  inc cx;
  cmp cx, R9w;
  jb @Loop;
  cmp bx, R10w;
  je @exit; // максимум в первой строке, обмен не нужен
  mov RCX, [R8+RBX*8];
  mov RAX, [R8+R10*8];
  mov [R8+RBX*8], RAX;
  mov [R8+R10*8], RCX;
@exit:
  xor dx, dx;
  fucomip ST(0), ST(1);
  mov cx, 1;
  fstp ST(0);
  cmovb ax, cx;
  cmovae ax, dx;
end;

procedure Make1RowPASM;
asm
// R8 - P, R9w - n, R10w - index
  mov R11, [R8 + R10*8];
  fld1;
  fld dword ptr [R11 + R10*4];
  mov BX, R10w;
  fdivp;
@Loop:
  fld st(0); // копируем d;
  inc BX;
  fmul dword ptr [R11 + RBX*4]; // в ST(0) - P[i]*d
  fstp dword ptr [R11 + RBX*4]; // созраняем ST(0) в P[i]
  cmp BX, R9w;
  jb @Loop;
  fstp st(0); // очищаем стек
end;

procedure MakeSingleNextRowsPASM;
asm
// R8 - P, R9w - n, R10w - index
  lea EBX, [R10d+1]; // BX - i
  xor RDX, RDX;
  cmp BX, R9w;
  jae @exit;
@Loop1:
  mov RCX, [R8 + RBX*8]; // RCX - S
  mov DX, R10w;
  fld dword ptr [RCX+R10*4];
@Loop2:
  fld st(0);
  inc DX;
  fmul dword ptr [R11+RDX*4];
  fsubr dword ptr [RCX+RDX*4];
  fstp dword ptr [RCX+RDX*4];
  cmp DX, R9W;
  jb @Loop2;
  fstp st(0);
  inc BX;
  cmp BX, R9w;
  jb @Loop1;
@exit:
end;

procedure ReverseSingleGJPASM;
asm
// R8 - P, R9w - n, R10w - index
  mov RDX, [R8+R10*8+8]; // RDX - R
  lea EBX, [R10d+1];
  fld dword ptr [RDX+R9*4];
@Loop:
  dec BX;
  mov R11, [R8+RBX*8];
  fld st(0);
  fmul dword ptr [R11+R10*4+4];
  fsubr dword ptr [R11+R9*4];
  fstp dword ptr [R11+R9*4];
  test BX, BX;
  jnz @Loop;
  fstp st(0);
end;