21.10.2020

Lazarus и AVX

Lazarus - это IDE а-ля Delphi, только под Free Pascal (FPC). Установил, на первый взгляд все было очень ничего. Внешний вид очень похож на Delphi 7.0 - самый удобный интерфейс для разработки. Тот интерфейс, что появился в более поздних версиях Delphi, лучше использовал площадь широкоформатного экрана, но был гораздо менее удобен.
Единственно, что удивило, в настройках интерфейса по умолчанию вся пунктуация (это называется "языковые символы") - красного цвета. Жесть.

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

С ассемблером оказалось сложнее. Споткнулся компилятор вот на таком фрагменте: 

procedure SolveTestSingleLESP_SSE(var P : TLESInfoP; LESNum : UInt64); // P - RCX, LESNum - RDX;
asm
// 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]; //
  ...

Параметры передаются в регистре RCX и RDX соответственно (это, оказывается, в Windows так принято, поэтому работает одинаково что в Delphi, что в FPC). И если с первыми строчками все нормально, то вот с последней компилятор забавно загнал. Он правильно распознал, что LESNum находится в RDX, но почему-то транслировал после замены переменной на регистр в mov R8, [R8*8+RDX].
Если переменную заменить явно на имя регистра, то транслирует нормально.
Второй затык случился на команде такого типа:

movss XMM0, dword ptr c_single1[0]; 

Почему-то компилятор выдал предупреждение, что он сгенерировал абсолютную 32-битную ссылку на константу, поэтому программа обязательно должна быть загружена в младшие четыре гигабайта ОЗУ!
Едрит-Мадрид! Это как и почему? Вроде же используется виртуальная память, которой все равно, где физически располагается переменная. И зачем специально генерировать абсолютную ссылку?
Естественно, в этом месте возникает ошибка. Наверное, это можно как-то порешать, используя какие-то специальные конструкции для ассемблера, имеющиеся в FPC. Но мне стало скучно и лень. Тем более, нормальной документации нету, что бы можно было легко докопаться до таких тонкостей.

Отдельно про debbuger. Я выяснил, что в Delphi он глючит, когда FPU (x87) занято данными более чем наполовину (просто данные верхней, самой новой половины стека не показывает). В Lazarus такого глюка нет. Просто потому, что там нет ресурсов FPU в отладчике как такового.
Ну и вообще там очень хиленькие и неудобные возможности по отладке в ассемблере по сравнению с Delphi. Регистров AVX там, кстати, тоже нет.
Синтаксис AT&T ассемблера отдельно достает. Нет, если впитал его с молоком матери, то нормально. Но у меня обратная ситуация. И если при разработке можно переключить ассемблер на синтаксис Intel, то отладчик всегда показывает только
AT&T синтаксис.
В целом Lazarus произвел впечатление полудоделанного и слегка подзаброшенного продукта.

В общем, плюнул я на это дело и решил попробовать запасной вариант: написать нужные мне процедуры на чистом ассемблере.
Идеальный для меня вариант - TASM. Но беда в том, что Embarcadero его давно забросил. И поэтому он поддерживает ровно тот же набор команд, что и компилятор Delphi: то есть AVX там нет.
Почитал в инете, вроде прямым наследником TASM, который более-менее развивается, является FASM. Это как бы и так, но не совсем. FASM настолько далеко ушел в своем развитии, что я не смог сгенерировать объектный файл с самыми простыми процедурами. Отказывается последняя версия FASM делать такой объектный файл. Хотя вроде долго смотрел в примеры, писал все по образцу, но не получается. Правда, не догадался попробовать сами примеры скачать и скомилировать... Надо будет как-нибудь ради прикола сделать.
Короче, заменил процедуры метками, благо, у меня все передается в регистрах, а локальные переменные мне не нужны. Видимо, если бы были нужны, надо было стек вручную создавать? Ну да бог с ним.

Быстренько скопировал код SSE, адаптировал его под AVX и на типе single (float), удивительное дело, все правильно сработало с первого раза, без отладки. Но, ради интереса, решил посмотреть под отладчиком. И не смог: при попытке зайти в процедуру с кодом под AVX отладчик Delphi просто "умирает".
Так что отладки под AVX у меня нету. С кодом под double переход был не таким гладким, в паре мест накосячил, пришлось повозится, выискиваю баги без отладчика. Но тем не менее все получилось.

Теперь осталось собрать тестер и оттестировать доступные мне процессоры. Но предварительно уже могу сказать: на PileDriver код на AVX немного медленнее, чем на SSE. Впрочем, это ожидаемо: давно известно, что PileDriver очень медленно выполняется обращение к памяти по 256 бит. А большинство команд AVX именно такие.
И даже команды FMA PileDriver не помогли. Работает ровно с той же скоростью, что и просто AVX. Но точность чуть выше, чем отдельно умножение + сложение. Впрочем, для метода Гаусса-Жордана повышение точности не очень существенное.
Наверное, если использовать команды AVX на половинных регистрах (XMM
вместо YMM), то может быть AVX слегка и обгонит SSE на моем FX-4350, особенно если FMA использовать. Но, мне кажется, это будет не совсем честно по отношению к другим процессорам 😉.

12.10.2020

Маленький размер

Заболел. В перерывах постельного режима решил ускорить работу СЛАУ как раз самого маленького размера: в 3 уравнения. Все же похожие системы очень часто встречаются в реальности. Практически все трехмерные расчеты так или иначе связаны примерно с этим размером. В то время как решение очень больших СЛАУ - это скорее очень редкие исключения. По крайней мере мне так кажется.
 
Начал, естественно, с FPU (x87). У него есть 8 вещественных регистров. А коэффициентов в СЛАУ 12. Впрочем, точно хранить нужно только элементы над верхней диагональю и свободные члены, которых всего 7. Но еще нужны свободные регистры для временных операций...
Короче, программирование матричных операций на регистрах, завязанных в стек ‒ та еще радость! И самое главное, не очень ясно, велик ли выигрыш от хранения коэффициентов в регистрах.

Попутно полдня мучил свой воспаленный мозг, пытаясь понять, что же не так то!? Через несколько часов выяснил, что "не так" - это глюк в дебаггере Delphi, который по крайней мере в 64-битном режиме неверно отображает заполненный более чем на половину регистровый стек FPU!
Да, видимо придется с Delphi все же уходить, как не жаль. Тем более, что Lazarus (точнее, Free Pascal) умеет вроде даже в AVX.

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

04.10.2020

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;