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;