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;

Про указатели

Итак, накидал тут намеднясь первую версию теста с заменой указателей на индексы в массиве. Что бы избавиться от длинных операций с памятью, заменив их вычислением адреса. Надо сказать, что 16 РОН тут решают, без них вся затея не стоила бы выделки.
В итоге получилось, что одна операция чтения из памяти 8 байт заменяется на чтение 2 байт, одно 32-битное умножение и сложение.

И вот здесь возникло у меня сомнение. А точно ли второе быстрее первого? Решил проверить. Чисто на паскале это вообще не так, с указателями работает гораздо быстрее. А вот на ассемблере все не так однозначно.
Код на указателях почти всегда чуть быстрее кода на индексах, малоуловимо, чуть выше погрешности измерения.

С памятью, конечно, ситуация не в пользу указателей. Если использовать указатели, то в худшем случае (СЛАУ из 3-х переменных типа float32) на хранение коэффициентов можно использовать 60% памяти, 40% уйдет на указатели. В случае же использования индексов можно под коэффициенты использовать не менее 83% памяти.
Так что даже не знаю, какой вариант выбрать. Как думаете?

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

@LOOP:
  inc BX;
  fmul dword ptr [R11 + RBX*4]; // вычисление адреса займет 2 такта
  ...
  jb @LOOP;

заменить таким
  lea RDX, [R11 + RBX*4];
@LOOP:
  add RDX, 4;// вычисление адреса займет 1 такт
  inc BX;
  fmul dword ptr [RDX];
...
  jb @LOOP;

Идея состояла в том, что бы немного ускорить перемещение от элемента к элементу. Так вот, такой вариант оказался  медленнее, по крайней мере для операций с плавающей точкой на моем процессоре. Честно говоря, несколько удивился этому.

64-бита, указатели и массивы

Итак, потихоньку переезжаю на 64-битный режим. Его преимущество очевидно - существенно больший, практически не ограниченный в обозримом будущем объем доступной оперативной памяти. Но есть и недостаток - размер указателей стал в 2 раза больше. 😀

Чем же это плохо? Плохо это становится, когда требуется интенсивная работа с указателями. Например, в случае двоичного дерева требуется на один узел хранить 2 указателя. Пусть в каждом узле мы храним 4-байтное целое.
Тогда на 32-битных системах полезная нагрузка будет всего лишь 33%, а вся остальная память уйдет на организацию структуры данных. При переходе на 64-разрядную систему полезной нагрузки будет всего лишь 16%, а основная часть памяти уйдет на структуру.

Казалось бы, какая разница? Ведь на 64-битных системах памяти, как говорится, хоть жопой жуй. Но это чисто теоретически. Практически же объем памяти в реальных системах ограничен как экономическими соображениями, так и техническими возможностями материнской платы. Мне кажется, в типичных настольных системах в настоящее время 16 ГиБ не часто встречается, а 32 ГиБ - очень редко. Поэтому если вам не хватало памяти в 32-битном режиме, не факт, что при переходе в 64-битный ее в реальности будет так много, что вообще не нужно будет думать об экономии.
Но основная причина все же не нехватка памяти из-за размера указателей, а снижение скорости работы, так как по шине данных нужно передавать в 2 раза больше данных, плюс снижается эффективность работы кэша. Допустим, у нас была 32-битная система с частотой памяти 1600 МГц, ее портировали на 64-битную и запустили на новой системе с частотой памяти 3200 МГц. И в результате (если в системе используется интенсивная работа с указателями), получили такую же производительность, как и на старой системе!

Поэтому, если программе требуется много работать с указателями и она умещается в 32-битное адресное пространство, то лучше ее в 32-битном режиме и оставить. Тем более, что современные 64-битные ОС позволяют без потери эффективности выполнять и 32-битные программы. Кстати, как не печально, в Windows прикладной программе по умолчанию доступно всего лишь половина из всего адресного пространства. Смешное ограничение, из разряда 640 Кбайт хватит всем.
Но если адресного пространства 32-битного режима не хватает и требуется максимальная эффективность по быстродействию и памяти, то можно попробовать заменить указатели на индексы массива, которые зачастую укладываются в меньшей размер, чем указатель. Ограничение такого подхода состоит в том, что трудно реализовать сильно меняющиеся в процессе работы программы структуры.

Ладно, вернемся к нашим баранам, то есть к тестированию производительности. Казалось бы, где там указатели, когда СЛАУ задается матрицей? Давайте посмотрим.
Набор матриц для СЛАУ можно описать разными способами. Самый простой:

TMatrices = array[1..k] of array[1..n] of array[1..n+1] of single;

Здесь для понятности массив индексируется с 1, хотя с 0 эффективнее. Впрочем, возможно, компилятор оптимизирует все это дело, но не факт. Здесь k - количество систем, n - размерность СЛАУ. Последний индекс на единицу больше, так как свободные члены СЛАУ хранятся просто в последнем столбце.
Кроме простоты, ничего хорошего в данном варианте нет. Главный недостаток - количество систем и их размер нужно знать заранее, до компиляции программы. От этого недостатка легко избавиться, просто перейдя к одномерному массиву:

TMatrices = array[0..Size-1] of single;

Создав указатель на этот тип, можно выделять в процессе работы столько памяти, сколько необходимо. Но в этом случае придется вручную рассчитывать индекс нужной ячейки. То есть если для первого варианта достаточно просто указать индексы q-ой системы, i-ой строки и j-го столбца, вот так

m[q, i, j]

а дальше компилятор сам сделает всю работы, то во втором случае придется чуть усложнить:

m[q*n*(n+1) + i*(n+1) + j]

С точки зрения вычислительной сложности оба варианта эквивалентны, по крайней мере для оптимизирующего компилятора. Во втором случае у нас появляется указатель, так как память выделяется динамически по мере необходимости, но он всего лишь один и никакой интенсивной работы с ним нет.
Но у обоих вариантов есть еще один недостаток: при решение СДАУ методом Гаусса-Жордана для снижения ошибки решения необходимо выделять главный элемент, что требует обмена строк, а это,  в общем-то, является очень медленной операцией.
Для решения проблемы можно было бы воспользоваться динамическими массивами, тогда объявление типа выглядело бы так:

TMatrices = array of array of array of single;

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

TRow = array[0..65535] of single;
PRow = ^TRow;
TMatrix = array[0..65535] of PRow;
PMatrix = ^TMatrix;
TMatrices = array[0..65535] of PMatrix;
PMatrices = ^TMatrices;

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

Поэтому попробую от указателей избавиться. Видится мне это пока примерно так:

TLESData = array[0..65535] of single;
PLESData = ^TLESData;
TLESRows = array[0..65535] of cardinal;
PLESRows = ^TLESRows;
TLESSet = array[0..65535] of PLESData;
PLESSet = ^TLESSets;

Здесь TLESData ‒ одномерный массив, как во втором, рассмотренном выше, варианте. Для того, что бы можно было быстро менять строки местами, используется тип TLESRows, в котором хранятся индексы первых элементов строк в массиве TLESData.
Так как в TLESRows хранятся 32-битные беззнаковые целые, он позволяет адресовать до 4 миллиардов (2³²) чисел в TLESData. Это позволит использовать до 16 ГиБ памяти под коэффициенты СЛАУ в случае FPU32, для других вещественных типов пропорционально больше. В случае СЛАУ из 3-х переменных это позволит хранить почти 358 млн. систем!
Адресация элементов матрицы, правда, в этом случае будет немного другой, чем во втором случае:

m[Rows[q*n + i] + j]

Мне кажется, такой вариант очень эффективно оптимизируется и будет быстрее, чем обращение просто по указателю.
Зачем же нужен третий тип, TLESSet? Дело в том, что производительность современных процессоров такова, что вполне вероятно топовые процессоры имеют или в ближайшем будущем будет иметь такую производительность, что смогут решить 358 млн. систем быстрее, чем за секунду! И уж тем более могут быть быстрее, когда речь идет про все ядра. И в таких случаях мне и потребуется несколько одномерных массивов по 16 ГиБ каждый в пределе, что бы можно было оттестировать такие монструозные процессоры.

Теперь остается попытаться реализовать все это на 64-битной версии Delphi. Прямо интересно, удастся ли это мне сделать или нет. В частности, на текущий момент та версия, что стоит у меня, даже в 64-битном режиме не позволяет создавать и обращаться к массивам размером больше 2 ГиБ. Естественно, на ассемблере эти ограничения мне нипочем, но не хотелось бы вообще все писать на нем. Впрочем, у меня есть идеи, как обойти это ограничение.
Хорошо хоть вроде бы память под динамические структуры выделяется полноценно. Но посмотрим.

Тестер. Первые результаты

Особенности работы протестированных процессоров в параллельном режиме я описал в предыдущем посте. Теперь посмотрим, на что способно одно ядро каждого процессора.

FPU32. СЛАУ из 3-х уравнений AMD FX-4350 решает 8.5 млн./сек. Это несколько отличается от значения, которое я приводил в самой первой статье. Объясняется просто: там я оценивал производительность, решая много раз одну и ту же СЛАУ, поэтому данные попали в кэш первого уровня и производительность была больше. В конечно же варианте теста производительность оценивается при решении разных систем.

i3-3227U обеспечивает 6.7 млн./сек, а i7-6700HQ - 13.6 млн./сек.

FPU64. СЛАУ из 3-х уравнений FX-4350 решает также 8.5 млн./сек., i3-3227U - 6.7 млн./сек., а вот i7-6700HQ настолько быстр, что в 32-разрядном режиме памяти под все системы не хватило, что бы время решения было примерно 1 секунду. И для система из 6 уравнений тоже.

Системы же из 12 уравнений этот процессор решает 792 тыс./сек., FX-4350 - 454 тыс./сек, i3-3227U - 407 тыс./сек.

FPU80. СЛАУ из 3-х уравнений FX-4350 решает также 3.5 млн./сек., а (сюрприз!) i3-3227U - 4.5 млн./сек. То есть маленький, крохотный процессор для нетбуков от Intel с тактовой частотой 1900 МГц кроет как бык овцу AMDшный процессор с частотой 4200МГц!  i7-6700HQ опять же оказался настолько быстр, что памяти под такие СЛАУ не хватило.

А вот системы из 6-ти уравнений процессоры решили так: FX-4350 - 603 тыс./сек, i3-3227U - 1 млн./сек., i7-6700HQ - 2 млн/сек.

Какой вывод из всего этого можно сделать? Архитектура Piledriver от AMD очень сильно проигрывала решениям Intel, несмотря на более высокие тактовые частоты, по крайней мере при операциях с плавающей точкой.

А что касается научных вычислений с максимально возможной точностью, то AMD находилась где-то ниже плинтуса.

Правда, в последних поколениях процессоров AMD исправилась и сейчас отстает как раз Intel. Но пока протестировать последние процессоры от AMD мне не довелось.

Общая оценка. Методика оценки очень проста. Решается набор СЛАУ разных размерностей от 3 до максимальной, на решение которой требуется не более 1 секунды. Увеличение размерности идет по геометрической прогрессии со знаменателем 2.

Процессоры сравниваются только на тех размерностях, которые у них совпадают. Разные типы данных имеют разный приоритет. Для текущей оценки я взял приоритет 0.5 для FPU32, 0.4 для FPU64 и 0.1 для FPU80.

Однопоточная производительность имеет приоритет 0.6, параллельная - 0.4. Если за единицу брать производительность процессора FX-4350, то получается вот такая картинка.

FX-4350 имеет базовую частоту 4200 МГц, частоту шины памяти 1600, i3-3227U - 1900 и 1333 МГц, i7-6700HQ - 2600 и 2133 МГц соответственно

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

Что дальше? Дальше нужно решить пару проблем.

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

В принципе, при решении реальных СЛАУ такие ситуации нужно просто отслеживать и выдавать предупреждение в ходе решения. Но не в случае замера производительности!

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

Вторая проблема состоит в том, что современные процессоры очень-очень быстры. А для точного замера производительности нужно решить достаточное количество систем. Я в текущем варианте беру такое количество, которое решается примерно за 1 секунду. В таком режиме погрешность измерения составляет примерно 5% (хотя точно научными методами я ее не оценивал), поэтому уменьшать время не хотелось бы.

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

Но в 32-разрядном режиме для маленьких систем у меня не хватает памяти для их размещения! То есть 2 гигабайт памяти, выделяемой программам Windows в 32-разрядном режиме не достаточно, так как время решения всех систем, помещающихся в эту память, меньше 1 секунды! Ситуация еще более усложняется при работе в параллельном режиме, так как там каждому ядру необходимо решить свой отдельный набор СЛАУ.

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

Таким образом, крайне актуальной задачей является переход на 64-разрядную версию программы. Там есть свои нюансы, но об этом в следующий раз.

Доделки тестера

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

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

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

Близкое к теоретическому ускорение можно получить, решаю очень сложную в вычислительном плане задачу, которая занимает мало памяти, например, целиком помещаясь  в кэше L1 или меньше.

Применительно же к моему методу тестирование параллельный режим состоит в том, что каждый процессор решает свой независимый набор систем линейных алгебраических уравнений (СЛАУ).

Кооперативный же режим предполагает, что процессоры совместно решают одну задачу. Для рассматриваемого теста - это когда одну большую СЛАУ решают сразу несколько процессоров. В принципе, метод Гаусса-Жордана относительно легко укладывается в этот режим, но все же гораздо сложнее, чем параллельный.

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

Но могу предположить, что в этом режиме ускорение будет еще меньше, чем в параллельном, так как конфликты при обращении к памяти будут чаще. Кроме того, часть данных может быть кэширована в L1/L2 одного ядра, когда к нему обратится второе, что скорее всего приведет к дополнительным задержкам. И тут может сыграть различная система организаций кэша в процессорах Intel и AMD. Вроде как раньше в AMD кэш был строго эксклюзивным, в то время как в Intel наоборот, инклюзивным. Но это не точно. Тем не менее, разница в стратегии кэширования может приводить к существенной разнице и в производительности.

Давайте посмотрим, что показали параллельные тесты.

Для начала мой родной, сильно постаревший и сильно бюджетный AMD FX-4350, имеющий 4 ядра. Для FPU32 и СЛАУ из 48 переменных параллельный режим дает ускорение примерно в 2.9 раза. По мере роста размерности ускорение падает до 2.3 для 768 переменных.

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

Однако такое предположение опровергается тем, что для FPU64, имеющему почти такую же производительность, но в 2 раза большее количество пересылок по шине памяти, ускорение совпадает почти один в один.

Ситуация меняется кардинально при работе с FPU80. Ускорение для СЛАУ 12 переменных составляет 3.5 и по мере увеличения размерности задачи растет, достигая 3.7 для 384 переменных! Предполагаю, что операции с расширенной точностью требуют на одну команду больше и на процессорах AMD выполняются достаточно медленно, снижая таким образом частоту конфликтов при обращении к памяти.

Таким образом, можно считать что на процессорах AMD поколения Piledriver эффективность параллельной работы составляет примерно 75% или чуть меньше при работе FPU.

Теперь посмотрим на мобильные решения от Intel. Начнем со старенького  i3-3227U, это двухядерный процессор с HT, с базовой частотой 1900 МГц с ограниченным теплопакетом, использовался в нетбуках. Этот процессор вообще выпрыгивает из штанов, показывая ускорение  в 2.6 для FPU32. То есть ускорение существенно больше, чем физических ядер! Вот здесь хорошо видно пользу HT. Фактически можно сказать, что HT на двухядерном процессоре добавляет еще пол ядра.

Больше никаких нюансов на этом процессоре нет, для всех типов данных максимальный коэффициент лежит в районе 2.4-2.7, при увеличении размерности коэффициент ускорения монотонно падает, за исключением очень больших размеров СЛАУ. Впрочем, такие размерности процессор решает настолько медленно, что особо я их не тестировал.

Еще один мобильный процессор Intel i7-6700HQ, поновее, 4 ядра, HT, частота 2600 МГц. Выступает чуть хуже предыдущего, из штанов практически не выпрыгивает: на малых размерностях ускорение может быть чуть больше 4, что также больше количества физических ядер. По мере увеличения размерности ускорение монотонно падает, сильно снижаясь после примерно 384 переменных. Видимо, это связано с выходом размера СЛАУ за пределы кзша L2, который у этого процессора равен 1 МБ.

Таким образом, Intel за счет технологии HT имеют некоторое преимущество над AMD, позволяя им практически на 90-100% использоваться физические ядра.

Впрочем, это не говорит о том, что технология AMD (называется CMT) хуже аналогичной Intel (которая называется похоже - SMT). Просто у AMD в Piledriver один исполнительный блок на два ядра. В каждом блоке два кластера для целочисленных операций, и лишь один - для операций с плавающей точкой. У Intel ситуация более честная. Тут скорее удивительно, что не смотря на то, что каждому ядру досталось лишь полкластера для операций с плавающей точкой, AMD все же показывает 75% эффективность параллельной работы.

Развлекаюсь, как могу

В этот раз могу так: взять любимый язык Pascal и на нем что-нибудь наваять. Наваять захотелось какую-нибудь мерилку производительности процессора. Люблю, я, понимаешь, всякими числами всё измерять.

Понятно, что мерилок таких понаделано и без меня не один вагон. Но все же хочется чего-то своего. Думал, думал, как мерить, и не придумал ничего лучше, чем решением систем линейных уравнений методом Гаусса-Жордана. Не знаю, почему именно этот метод мне понравился, но вроде норм.

Такой подход хорош тем, что можно замерить не только систему команд SISD, но и SIMD, смотря как реализовывать. SIMD мне давно хотелось попробовать, но все руки никак не доходили, да и там на x86/x64 такой зоопарк из инструкций, что просто так, без бутылки, к ним лучше и не лезть.

Ну, то есть вы поняли уже, да? Я не только Pascal люблю, но и Assembler в лайтой версии тоже. Видимо, детская травма суровым советским бытом: калькулятор MK-61 и вот всё это вот 😉.

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

Дальше мне стало скучно и я переписал решение на ассемблере. Сначала начал для чисел расширенной точности, в терминах Pascal это тип extended или, для краткости, буду здесь его называть FP80 (потому что занимает 80 бит).

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

Думаю, для более коротких вещественных типов выигрыш будет, наверное, еще меньше. Но все же реализовал ради забавы. И оказалось, что я ошибался. Для типа с двойной точноностью (double или FP64) выигрыш оказался 1.76 раза для маленьких и 3.44 раза для больших систем (1536 неизвестных).

Для типа одинарной точности (single или FP32) выигрыш составил 1.7 раза для маленьких и 5 раз для больших систем (1536 неизвестных). Конечно, я город таким достижением 😀, хотя на самом деле это не моя заслуга, а Delphi, так как ходят слухи, что её последние версию не очень хороши в оптимизации.

Кстати, заметна разница в производительности между FP80 и остальными. Это связано с тем, что в FPU нет команд для операций с FP80 прямо из памяти. То есть приходится загружать оба операнда в регистры, проводить операцию между регистрами, потом сохранять. Примерно так:

  fld tbyte ptr [eax];
  fmul st(0), st(1);
  fstp tbyte ptr [eax]; 
А для FP32 это выглядит чуть по-другому:
  fld st(0);  
  fmul dword ptr [eax];
  fstp dword ptr [eax];Вроде команд столько же, но работают они гораздо быстрее, такая уж особенность архитектуры: при выполнении операции с операндом в памяти загрузка из памяти и расчет значения производятся, похоже, одновременно.
Ну и сказывается, конечно, то, что 10 байт плохо выравниваются, что снижает эффективность и кэширования, и операций с памятью.

Предварительные результаты: на стареньком бюджетном AMD FX-4350 частотой 4.2 ГГц и памятью на 1600 МГц процессор решает 3.8 млн. систем FP80 из 3-х линейных уравнений за 1 секунду, 10.7 млн систем FP64 и 10.8 млн. FP32.

Правда, сложность решения пропорциональна кубу количества переменных, так что максимум, что можно решить за одну секунду, составляет в районе несколько сотен переменных FP80 (< 700), и чуть больше 1000 для FP32.

Много это или мало? Смотря для чего. Скажу так, что 30 лет назад о таком и мечтать было нельзя, причем не только в г. Новокузнецке. Но, например, при расчетах, связанных с моделированием реального мира, это не много. Например, при расчете модели трехмерного мира модным ныне методом трассировки лучей, где используются не совсем такие, но близкие к этому методы, такая производительность даже близко не позволит подойти к получению высококачественной картинки в реальном времени.

Ну а я пока продолжу развлечение. Если будет что интересное - напишу.

Шифрация. Продолжение.

Ого! Как активно я пишу в блог-то. ))) С прошлого сообщения про шифрацию прошел ровно год.

За этот год мне в голову пришли некоторые улучшения в описанный ранее алгоритм.

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

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

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

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

В четвертых, это вариант использование ЛКГ. Считается, что младшие биты ЛКГ слишком легко предсказуемы и статически связаны между собой, поэтому лучше использовать для генерации случайных значений биты подальше от правого края числа.

Но меня смутило другое соображение. Пусть у нас имеется параметр a, который является множителем ЛКГ. С этим параметром и заданным ключом X, которое является начальным значением для ЛКГ, у нас получится ряд случайных значений R1,..., Rn, которое мы получили, взяв справа некоторое количество бит из сгенерированных значений ЛКГ X1, ..., Xn. Не важно, возьмем ли мы их совсем справа, или с каким-то смещением.

Теперь, допустим, мы захотели увеличить в два раза разрядность ЛКГ, для существенного повышения криптостойкость. Если мы это сделаем таким образом: добавив к a слева такое же по разрядности число b, то есть ba, а к ключу - Q, получив QX, и используя тот же алгоритм выделения случайных значений, мы получим тy же последовательность R1,..., Rn, что и ранее! То есть на самом деле криптостойкость не увеличилась, и можно раскрыть даже полным перебором с небольшими затратами времени более криптостойкий алгоритм. Возможно, использую даже меньшую разрядность, что исходные a и X.

Поэтому для шифрации надо брать самые старшие разряды из сгенерированного ЛКГ числа. Не знаю, может в этом случае тоже есть какие-то недостатки, я пока их не заметил.

Ну и наконец-то, не прошло и года, как я таки накидал первую, корявую версию на 64-битном ЛКГ генераторе. На всякий случай оставлю ее здесь, а то когда-то давно я уже делал что-подобное, менее совершенное. И все пропала из-за умершего HDD. И да, на текущий момент я уверен, что мой метод абсолютно надежен. Надеюсь, сами догадаетесь почему. ;)