15.12.2020 17:34
rebuilder 38 3200 Источник


Приветствую всех любителей алгоритмов. Хочу Вам поведать о своих изысканиях на тему сортировок в целом и углубиться в рассмотрение radix сортировки.

Будучи разработчиком с многолетним стажем, всё чаще стал сталкиваться со странной тенденцией в разработке программного обеспечения:

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

Полагаю это связано с общей идеей отдать предпочтение быстрому программированию с использованием всё более мощных фреймворков и сверхвысокоуровневых / скриптовых языков программирования. Языки подобные Ruby или Python невероятно удобны для разработчика. Множество 'синтаксического сахара', я бы даже сказал 'Мёда', ускоряют разработку в разы, если не на порядки, но какой ценой. Как пользователя, меня напрягает низкая тепловая эффективность кода, про объёмы потребляемой памяти просто промолчу, однако главный ресурс человечества – время. Оно бесследно исчезает в бесконечных абстракциях, похищается анализаторами кода, вычищается умными сборщиками мусора. Я не призываю вернуться в прошлое, отказавшись от благ современной разработки, писать “дорогой” код, лишь предлагаю задуматься над возможным устранением “бутылочных горлышек” производительности там, где это возможно в типовых задачах. Часто этого можно достичь, оптимизировав высоконагруженные участки кода.

Как одну из базовых задач оптимизации можно выделить сортировку. Тема настолько исследована, вдоль и поперёк, что, казалось бы, на этом пути трудно что-то обнаружить интересного. Однако мы попробуем.

Сортировать малые массивы (менее миллиона элементов) мы не станем. Даже, если это делать крайне неэффективно, почувствовать просадки достаточно сложно, так как они нивелируются производительностью современного оборудования. Другое дело большие объёмы данных (миллиарды элементов), от грамотного подбора алгоритма очень сильно варьируется скорость выполнения.

Все алгоритмы, основанные на сравнениях, в общем случае решают задачу сортировки не лучше чем O(n * Log n). При больших n эффективность быстро снижается и изменить эту ситуацию не представляется возможным. Эту тенденцию можно исправить, отказавшись от методов, основанных на сравнении элементов. Наиболее перспективным мне видится алгоритм Radix sort. Его вычислительная сложность O(k * n), где k количество проходов по массиву. Если n достаточно велико, а k наоборот очень мало, то данный алгоритм выигрывает у O(n * Log n).
В современной архитектуре процессоров k практически всегда сводится к числу байт сортируемого числа. Так например, для DWord (int) k = 4, что весьма не много. Это некоторая “потенциальная яма” вычислений, которая обусловлена несколькими факторами:

  • Регистры процессора заточены под 8-битные операторы на аппаратном уровне
  • Используемый буфер подсчёта в алгоритме как раз укладывается в одну линию кэша L1 –процессора. (256 * 4-байтных числа)

В истинности этого утверждения Вы можете попробовать убедиться самостоятельно. Однако на текущий момент времени делить по биту самый оптимальный вариант. Не исключаю, что когда кэш L1 процессоров разрастётся до 256 КБайт, более выгодным станет вариант делить по границе 2-байта.

Эффективная реализация сортировки это не только алгоритм, но и тонкая инженерная задача по оптимизации кода.

В данном решении алгоритм состоит из нескольких этапов:

  1. Выделение памяти под временный массив, на фоне общего времени исполнения не самая дешёвая операция, как вариант можно вынести за пределы функции и подавать параметром
  2. Эффективный подсчёт по всем битам сразу за один проход с максимальной эффективностью кодогенерации
  3. Расчёт смещений, также для всех проходов сразу
  4. Собственно пошаговая сортировка по битам от младших к старшим (LSD), в виде отдельной функции

Алгоритм LSD применяем как более быстрый (по крайней мере, в моём варианте) за счёт более ровной обработки при различных флуктуациях входных данных.
Исходный полностью отсортированный массив является наихудшим случаем для алгоритма, так как данные всё равно будут полностью сортироваться. Зато на случайных или перемешанных данных Radix sort чрезвычайно эффективен.

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

typedef struct TNode {
  //unsigned long long key;
  unsigned int key;
  //unsigned short key;
  //unsigned  char key;
  unsigned int value;
  //unsigned int value1;
  //unsigned int value2;
} TNode;

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

Раннее я уже писал заметку с реализацией Radix sort на паскаль, однако “сегрегация” языков программирования набирает невиданные ранее темпы среди завсегдатаев данного ресурса. Поэтому код для данной статьи я решил частично переписать на 'си' как более 'правильный' распространённый язык в сообществе программистов (Хотя на выходе ассемблерная кодогенерация с Pascal местами практически идентичная). При сборке рекомендую применять ключ –O2 или выше.

В отличие от прошлой реализации, применив указатели в адресации цикла вместо операций побитного сдвига, удалось получить ещё большую прибавку в 1-2% к скорости алгоритма.

C 
#include <stdio.h>
#include <omp.h>

#include <time.h>
#include <windows.h>
#include <algorithm>
//=============================================================
typedef struct TNode {
  //unsigned long long key;
  unsigned int key;
  //unsigned short key;
  //unsigned  char key;
  unsigned int value;
  //unsigned int value1;
  //unsigned int value2;
} TNode;
//=============================================================
void RSort_step(TNode *source, TNode *dest, unsigned int n, unsigned int *offset, unsigned char sortable_bit)
{
  unsigned char *b = (unsigned char*)&source[n].key + sortable_bit;
  TNode *v = &source[n];
  while (v >= source)
  {
    dest[--offset[*b]] = *v--;
    b -= sizeof(TNode);
  }
}
//=============================================================
void RSort_Node(TNode *m, unsigned int n)
{
  // Выделяем память под временный массив
  TNode *m_temp = (TNode*)malloc(sizeof(TNode) * n);

  // Заводим массив корзин
  unsigned int s[sizeof(m->key) * 256] = {0};
  // Заполняем массив корзин для всех разрядов
  unsigned char *b = (unsigned char*)&m[n-1].key;
  while (b >= (unsigned char*)&m[0].key)
  {
    for (unsigned int digit=0; digit< sizeof(m->key); digit++)
    {
      s[*(b+digit)+256*digit]++;
    }
    b -= sizeof(TNode);
  }

  // Пересчитываем смещения для корзин
  for (unsigned int i = 1; i < 256; i++)
  {
    for (unsigned int digit=0; digit< sizeof(m->key); digit++)
    {
      s[i+256*digit] += s[i-1+256*digit];
    }
  }

  // Вызов сортировки по битам от младших к старшим (LSD)
  for (unsigned int digit=0; digit< sizeof(m->key); digit++)
  {
    RSort_step(m, m_temp, n-1, &s[256*digit] ,digit);
    TNode *temp = m;
    m = m_temp;
    m_temp = temp;
  }

  // Если ключ структуры однобайтовый, копируем отсортированное в исходный массив
  if (sizeof(m->key)==1)
  {
    TNode *temp = m;
    m = m_temp;
    m_temp = temp;
    memcpy(m, m_temp, n * sizeof(TNode));
  }

  free(m_temp);
}
//=============================================================
int main()
{
  unsigned int n=10000000;

  LARGE_INTEGER frequency;
  LARGE_INTEGER t1, t2, t3, t4;
  double elapsedTime;

  TNode *m1 = (TNode*)malloc(sizeof(TNode) * n);
  TNode *m2 = (TNode*)malloc(sizeof(TNode) * n);
  srand(time(NULL));

  for (unsigned int i=0; i<n; i++)
  {
      m1[i].key = rand()*RAND_MAX+rand();
      m2[i].key = m1[i].key;
  }

  QueryPerformanceFrequency(&frequency);
  QueryPerformanceCounter(&t1);

  RSort_Node(m1, n);

  QueryPerformanceCounter(&t2);
  elapsedTime=(float)(t2.QuadPart-t1.QuadPart)/frequency.QuadPart;
  printf("The RSort:          %.5f seconds\n", elapsedTime);

  QueryPerformanceFrequency(&frequency);
  QueryPerformanceCounter(&t3);

  std::sort(m2, m2+n,[](const TNode &a, const TNode &b){return a.key < b.key;});

  QueryPerformanceCounter(&t4);
  elapsedTime=(float)(t4.QuadPart-t3.QuadPart)/frequency.QuadPart;
  printf("The std::sort:      %.5f seconds\n", elapsedTime);

  for (unsigned int i=0; i<n; i++) 
  {
    if (m1[i].key!=m2[i].key) 
    {
      printf("\n\n!!!!!\n");
      break;
    }
  }

  free(m1);
  free(m2);
  return 0;
}


Pascal 
program SORT;
uses
  SysUtils, Windows;
//=============================================================
type TNode = record
     key : Longword;
     //value : Longword;
end;
type ATNode = array of TNode;
//=============================================================
procedure RSort_Node(var m: array of TNode);
//------------------------------------------------------------------------------
    procedure Sort_step(var source, dest: array of TNode; len : Longword; offset: PLongword; const num: Byte);
        var b : ^Byte;
            v : ^TNode;
        begin
          b:=@source[len];
          v:=@source[len];
          inc(b,num);
          while v >= @source do
          begin
            dec(offset[b^]);
            dest[offset[b^]] := v^;
            dec(b,SizeOf(TNode));
            dec(v);
          end;
        end;
//------------------------------------------------------------------------------
var // Объявляем массив корзин первым, для выравнивания на стеке
    s: array[0..1023] of Longword =(
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
    i : Longword;
    b : ^Byte;
    p : Pointer;
begin

  GetMem(p, SizeOf(TNode)*Length(m));

  // Подсчёт корзин
  b:=@m[High(m)];
  while (b >= @m[0]) do
  begin
    Inc(s[(b+3)^+256*3]);
    Inc(s[(b+2)^+256*2]);
    Inc(s[(b+1)^+256*1]);
    Inc(s[(b+0)^+256*0]);
    dec(b,SizeOf(TNode));
  end;

  // Пересчёт смещений для корзин
  for i := 1 to 255 do                             
  begin
    Inc(s[i+256*0], s[i-1+256*0]);
    Inc(s[i+256*1], s[i-1+256*1]);
    Inc(s[i+256*2], s[i-1+256*2]);
    Inc(s[i+256*3], s[i-1+256*3]);
  end;

  // Вызов сортировки по битам от младших к старшим
  Sort_step(m, ATNode(p), High(m), @s[256*0], 0);                  
  Sort_step(ATNode(p), m, High(m), @s[256*1], 1);
  Sort_step(m, ATNode(p), High(m), @s[256*2], 2);
  Sort_step(ATNode(p), m, High(m), @s[256*3], 3);

  FreeMem(p);
end;
//=============================================================
 procedure test();
  const
    n = 10000000;
  var
    m1: array of TNode;  
    i,j,k,j1,j2: Longword;

    iCounterPerSec: TLargeInteger;
    T1, T2: TLargeInteger; //значение счётчика ДО и ПОСЛЕ операции  
  begin
    SetLength(m1,n);

    for i := 0 to n - 1 do
    begin
      m1[i].key := Random(65536 * 65536);
    end;  

  QueryPerformanceFrequency(iCounterPerSec);//определили частоту счётчика
  QueryPerformanceCounter(T11); //засекли время начала операции

  RSort_Node(m1);
    
  QueryPerformanceCounter(T22);//засекли время окончания
  WRITELN('1='+FormatFloat('0.0000', (T22 - T11)/iCounterPerSec) + ' sec.');//вывели количество секунд на выполнение операции

    SetLength(m, 0);
  end;
  //------------------------------------------------------------------------------
begin
  test();
  Readln();
  exit;
end.


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

А если хочется ещё немного быстрее?

Дальнейшим логичным шагом, как мне виделось, использовать видеокарту. На просторах интернета я встречал множество рассуждений о том, что Radix sort отлично параллелился на множестве ядер видеокарты (чуть ли не лучший способ сортировки не видеокарте).

Поддавшись искушению получить дополнительную производительность, реализовал несколько вариантов сортировки на знакомом мне OpenCL. К сожалению у меня нет возможность проверить реализацию на топовых видеокартах, но на моей GEFORCE GTX 750 TI алгоритм проигрывал однопоточной реализации на CPU, за счёт того, что данные нужно отправить на видеокарту, а потом забрать обратно. Если не гонять данные по шине туда-сюда скорость была бы приемлемой, но всё равно не ‘фонтан’. Есть ещё одно замечание. В OpenCL нити выполнения не синхронны (рабочие группы выполняются в произвольном порядке, насколько мне известно, в CUDA это не так, поправьте, кто знает), что мешает написать в данном случае более эффективный код.

Код из серии: можно ли создать троллейбус… обработку в OpenCL на Дельфи... но зачем?
Переписывать на 'Си' поленился, выкладываю как есть.
program project1;

uses Cl, SysUtils, Windows, Math;
//------------------------------------------------------------------------------
function OCL_Get_Prog(context: cl_context; pdevice_id: Pcl_device_id; name: PChar; S:AnsiString):cl_kernel;
var   Tex:PCHAR;
      Len:QWORD;
      PLen:PQWORD;
      Prog:cl_program;
      kernel:cl_kernel;
      Ret:cl_int;
begin
   Tex:=@S[1];
   Len:=Length(S);
   PLen:=@LEN;
   prog:=nil;
   kernel:=nil;

     // создать бинарник из кода программы
     prog:= clCreateProgramWithSource(context, 1, @Tex, @Len, ret);
     if CL_SUCCESS<>ret then writeln('clCreateProgramWithSource Error ',ret);
     // скомпилировать программу
     ret:= clBuildProgram(prog, 1, pdevice_id, nil, nil, nil);
     if CL_SUCCESS<>ret then writeln('clBuildProgram            Error ',ret);
     // создать объект ядра для исполнения программы
     kernel:= clCreateKernel(prog, name, ret);
     if  ret<>CL_SUCCESS then writeln('clCreateKernel            Error ',ret);
     // удаляем программу
     clReleaseProgram(prog);
     OCL_Get_Prog:=kernel;
end;
//------------------------------------------------------------------------------
var
   context:cl_context;
   kernel1, kernel2, kernel3, kernel4 :cl_kernel;
   Ret:cl_int;

   valueSize : QWord =0;
   s0 : PChar;

   platform_id:cl_platform_id;
   ret_num_platforms:cl_uint;
   ret_num_devices:cl_uint;
   device_id:cl_device_id;
   command_queue:cl_command_queue;
   S1,S2,S3,S4 : AnsiString;
   memobj1, memobj2, memobj3 :cl_mem;
   mem:Array of LongWord;
   size:LongWord;
   g_works, l_works :LongInt;

   iCounterPerSec: TLargeInteger;
   T1, T2, T3, T4: TLargeInteger; //значение счётчика ДО и ПОСЛЕ операции
   i,j,step:LongWord;

procedure exchange_memobj(var a,b:cl_mem);
var c:cl_mem;
begin
  c:=a;
  a:=b;
  b:=c;
end;

begin
   //---------------------------------------------------------------------------
   S1 :=
   // Шаг сортировки 1 (Oбнуление массива счётчика)
   '__kernel void sort1(__global uint* m1) {'+
   ' uint g_id = get_global_id(0);'+
   ' m1[g_id] = 0;'+
   '}';
   //---------------------------------------------------------------------------
   S2 :=
   // Шаг сортировки 2 (Подсчёт индексов по блокам)
   '__kernel void sort2(__global uint* m1, __global uint* m2, const uint len, const uint digit) {'+
   ' uint g_id = get_global_id(0);'+
   ' uint size = get_global_size(0);'+
   ' uint a = g_id / len;'+
   ' uchar key = (m1[g_id] >> digit);'+
   ' atomic_inc(&m2[key]);'+
   ' atomic_inc(&m2[256 * a + key + 256]);'+
   '}';
   //---------------------------------------------------------------------------
   S3 :=
   // Шаг сортировки 3 (Расчёт смещения)
   '__kernel void sort3(__global uint* m1, const uint len) {'+
   ' uint l_id = get_global_id(0);'+

   ' for (uint i = 0; i < 8; i++) {'+
   '  uint offset = 1 << i;'+
   '  uint add = (l_id>=offset) ? m1[l_id - offset] : 0;'+
   '  barrier(CLK_GLOBAL_MEM_FENCE);'+
   '  m1[l_id] += add;'+
   '  barrier(CLK_GLOBAL_MEM_FENCE);'+
   ' }'+

   ' for (uint i = 1; i < 1024; i++) {'+
   '  m1[i*256+l_id + 256] += m1[(i-1)*256+l_id + 256];'+
   ' }'+

   '  barrier(CLK_GLOBAL_MEM_FENCE);'+

   ' if (l_id>0) {'+
   '  for (uint i = 0; i < 1024; i++) {'+
   '   m1[i*256+l_id + 256] += m1[l_id-1];'+
   '  }'+
   ' }'+

   '}';

   //---------------------------------------------------------------------------
   S4 :=
   // Шаг сортировки 4
   '__kernel void sort4(__global uint* m1, __global uint* m2, __global uint* m3, const uint digit, const uint len) {'+
   ' uint g_id = get_global_id(0);'+
   ' for (int i = len-1; i >= 0; i--) {'+      // обратить внимание цикл обратный!
   '  uchar key = (m1[g_id*len+i] >> digit);'+
   '  m2[--m3[g_id*256 + key + 256]] = m1[g_id*len+i];'+
   ' }'+
   '}';
   //---------------------------------------------------------------------------

   // Получаем доступные платформы
   ret := clGetPlatformIDs(1,@platform_id,@ret_num_platforms);
   if CL_SUCCESS<>ret then writeln('clGetPlatformIDs          Error ',ret);
   // Получаем доступные устройства
   ret := clGetDeviceIDs(platform_id, CL_DEVICE_TYPE_GPU, 1, @device_id, @ret_num_devices);
   if CL_SUCCESS<>ret then writeln('clGetDeviceIDs            Error ',ret);

   clGetDeviceInfo(device_id, CL_DEVICE_NAME, 0, nil, valueSize);
   GetMem(s0, valueSize);
   clGetDeviceInfo(device_id, CL_DEVICE_NAME, valueSize, s0, valueSize);
   Writeln('DEVICE_NAME:  '+s0);
   FreeMem(s0);

   // Создаём контекст для устройства
   context:= clCreateContext(nil, 1, @device_id, nil, nil, ret);
   if CL_SUCCESS<>ret then writeln('clCreateContext           Error ',ret);

   // Создаём очередь команд для контекста и устройства
   command_queue := clCreateCommandQueue(context, device_id, 0, ret);
   if CL_SUCCESS<>ret then writeln('clCreateContext           Error ',ret);

   //-------------------------------------------------------------
   kernel1 := OCL_Get_Prog(context, @device_id, 'sort1', S1);
   kernel2 := OCL_Get_Prog(context, @device_id, 'sort2', S2);
   kernel3 := OCL_Get_Prog(context, @device_id, 'sort3', S3);
   kernel4 := OCL_Get_Prog(context, @device_id, 'sort4', S4);
   //-------------------------------------------------------------

   size:=256*256*16*10;

   g_works := size;
   l_works := 256;

   Randomize;
   SetLength(mem, size);

   for i:=0 to size-1 do mem[i]:=random(256*256*256*256);

   // Создаём буферы для ввода вывода данных
   memobj1 := clCreateBuffer(context, CL_MEM_READ_WRITE, size * sizeof(cl_uint), nil, ret);
   if  ret<>CL_SUCCESS then writeln('clCreateBuffer1            Error ',ret);
   memobj2 := clCreateBuffer(context, CL_MEM_READ_WRITE, size * sizeof(cl_uint), nil, ret);
   if  ret<>CL_SUCCESS then writeln('clCreateBuffer2            Error ',ret);
   memobj3 := clCreateBuffer(context, CL_MEM_READ_WRITE, (256+256*1024) * sizeof(cl_uint), nil, ret);
   if  ret<>CL_SUCCESS then writeln('clCreateBuffer3            Error ',ret);

   QueryPerformanceFrequency(iCounterPerSec); //определили частоту счётчика
   QueryPerformanceCounter(T1); //засекли время начала операции

   // Записываем данные в буфер
   ret := clEnqueueWriteBuffer(command_queue, memobj1, CL_TRUE, 0, size * sizeof(cl_int), @mem[0], 0, nil, nil);
   if  ret<>CL_SUCCESS then writeln('clEnqueueWriteBuffer      Error ',ret);

   QueryPerformanceCounter(T2);//засекли время окончания
   Writeln('write       '+FormatFloat('0.0000', (T2 - T1)/iCounterPerSec) + ' second.');//вывели количество секунд на выполнение операции

   QueryPerformanceFrequency(iCounterPerSec); //определили частоту счётчика
   QueryPerformanceCounter(T3); //засекли время начала операции


   for step:=0 to 3 do
   begin

   //-------------------------------------------------------------
   QueryPerformanceFrequency(iCounterPerSec); //определили частоту счётчика
   QueryPerformanceCounter(T1); //засекли время начала операции
   //-------------------------------------------------------------
   // 1 ШАГ (очитска массива)
   ret:= clSetKernelArg(kernel1, 0, sizeof(cl_mem), @memobj3 );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_1_1            Error ',ret);

   i := 256+256*1024;
   ret:= clEnqueueNDRangeKernel(command_queue, kernel1, 1, nil, @i, nil, 0, nil, nil);
   if  ret<>CL_SUCCESS then writeln('clEnqueueNDRangeKernel_1    Error ',ret);
   //-------------------------------------------------------------

   clFinish(command_queue); // Ожидаем завершения всех команд в очереди
   QueryPerformanceCounter(T2);  //засекли время окончания
   Writeln('step 1      '+FormatFloat('0.0000', (T2 - T1)/iCounterPerSec) + ' second.');//вывели количество секунд на выполнение операции

   QueryPerformanceFrequency(iCounterPerSec); //определили частоту счётчика
   QueryPerformanceCounter(T1);  //засекли время начала операции
   //-------------------------------------------------------------
   // 2 ШАГ (заполнение счётчиков)

   ret:= clSetKernelArg(kernel2, 0, sizeof(cl_mem), @memobj1 );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_2_1            Error ',ret);
   ret:= clSetKernelArg(kernel2, 1, sizeof(cl_mem), @memobj3 );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_2_2            Error ',ret);
   j := size div (1024);
   ret:= clSetKernelArg(kernel2, 2, sizeof(j), @j );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_2_3            Error ',ret);

   j:=step*8;
   ret:= clSetKernelArg(kernel2, 3, sizeof(j), @j );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_2_4            Error ',ret);

   i := size;
   ret:= clEnqueueNDRangeKernel(command_queue, kernel2, 1, nil, @i, nil, 0, nil, nil);
   if  ret<>CL_SUCCESS then writeln('clEnqueueNDRangeKernel_2    Error ',ret);

   //-------------------------------------------------------------
   clFinish(command_queue); // Ожидаем завершения всех команд в очереди
   QueryPerformanceCounter(T2);//засекли время окончания
   Writeln('step 2      '+FormatFloat('0.0000', (T2 - T1)/iCounterPerSec) + ' second.');//вывели количество секунд на выполнение операции
   QueryPerformanceFrequency(iCounterPerSec);//определили частоту счётчика
   QueryPerformanceCounter(T1); //засекли время начала операции
   //-------------------------------------------------------------
   // 3 ШАГ (пересчёт смещений)

   ret:= clSetKernelArg(kernel3, 0, sizeof(cl_mem), @memobj3 );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_3_1            Error ',ret);

   j := size;
   ret:= clSetKernelArg(kernel3, 1, sizeof(j), @j );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_3_3            Error ',ret);

   j := 256;
   ret:= clEnqueueNDRangeKernel(command_queue, kernel3, 1, nil, @j, @j, 0, nil, nil);
   if  ret<>CL_SUCCESS then writeln('clEnqueueNDRangeKernel_3    Error ',ret);

   //-------------------------------------------------------------
   clFinish(command_queue); // Ожидаем завершения всех команд в очереди
   QueryPerformanceCounter(T2);//засекли время окончания
   Writeln('step 3      '+FormatFloat('0.0000', (T2 - T1)/iCounterPerSec) + ' second.');//вывели количество секунд на выполнение операции
   QueryPerformanceFrequency(iCounterPerSec);//определили частоту счётчика
   QueryPerformanceCounter(T1); //засекли время начала операции
   //-------------------------------------------------------------

   // 4 ШАГ (сортировка)
   ret:= clSetKernelArg(kernel4, 0, sizeof(cl_mem), @memobj1 );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_4_1            Error ',ret);
   ret:= clSetKernelArg(kernel4, 1, sizeof(cl_mem), @memobj2 );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_4_2            Error ',ret);
   ret:= clSetKernelArg(kernel4, 2, sizeof(cl_mem), @memobj3 );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_4_3            Error ',ret);

   j:=step*8;
   ret:= clSetKernelArg(kernel4, 3, sizeof(j), @j );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_4_4            Error ',ret);

   j := size div (1024);
   ret:= clSetKernelArg(kernel4, 4, sizeof(j), @j );
   if  ret<>CL_SUCCESS then writeln('clSetKernelArg_4_5            Error ',ret);

   i := (1024);  // глобально ядер
   ret:= clEnqueueNDRangeKernel(command_queue, kernel4, 1, nil, @i, nil, 0, nil, nil);
   if  ret<>CL_SUCCESS then writeln('clEnqueueNDRangeKernel_4    Error ',ret);
   clFinish(command_queue);
   // Меняем местами указатели на массивы после каждого шага
   // Результат останется в memobj1
   exchange_memobj(memobj2, memobj1);
   //-------------------------------------------------------------
   clFinish(command_queue); // Ожидаем завершения всех команд в очереди
   QueryPerformanceCounter(T2);//засекли время окончания
   Writeln('step 4      '+FormatFloat('0.0000', (T2 - T1)/iCounterPerSec) + ' second.');//вывели количество секунд на выполнение операции
   //-------------------------------------------------------------
   end;

   QueryPerformanceCounter(T4);//засекли время окончания
   Writeln('all not R/W '+FormatFloat('0.0000', (T4 - T3)/iCounterPerSec) + ' second.');//вывели количество секунд на выполнение операции


   QueryPerformanceFrequency(iCounterPerSec);//определили частоту счётчика
   QueryPerformanceCounter(T1); //засекли время начала операции

   // Считываем данные из буфера
   ret:= clEnqueueReadBuffer(command_queue, memobj1, CL_TRUE, 0, size * sizeof(cl_int), @mem[0], 0, nil, nil);
   if  ret<>CL_SUCCESS then writeln('clEnqueueReadBuffer       Error ',ret);

   QueryPerformanceCounter(T2);//засекли время окончания
   Writeln('Read        '+FormatFloat('0.0000', (T2 - T1)/iCounterPerSec) + ' second.');//вывели количество секунд на выполнение операции


    // Освобождаем ресурсы
   clReleaseMemObject(memobj1);           // Высвобождаем память буфера
   clReleaseMemObject(memobj2);
   clReleaseMemObject(memobj3);
   clReleaseKernel(kernel1);              // Высвобождаем объект исполнения
   clReleaseKernel(kernel2);
   clReleaseKernel(kernel3);
   clReleaseKernel(kernel4);
   clReleaseCommandQueue(command_queue);  // Высвобождаем очередь
   clReleaseContext(context);             // Высвобождаем контекст устройства
   //-------------------------------------------------------------
   SetLength(mem, 0);
   readln;
end.


Остался ещё один не опробованный вариант – многопоточность. Уже используя только голый C и библиотеку OpenMP, я решил узнать какой эффект даст использование множества ядер CPU.

Изначально идея заключалась в разбивке исходного массива на равные части, передача их в отдельные потоки, а затем склейка их слиянием (merge sort). Сортировка шла не плохо, а вот слияние сильно замедляло всю конструкцию, каждая склейка равносильна дополнительному проходу по массиву. Эффект был хуже чем работа в один поток. Забраковал реализацию, как не практичную.

В итоге применил параллельную сортировку очень похожую на ту, что применялась на GPU. Вот с ней уже всё получилось на много лучше.

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

C 
#include <stdio.h>
#include <omp.h>
//=============================================================
typedef struct TNode {
  //unsigned long long key;
  unsigned int key;
  //unsigned short key;
  //unsigned  char key;
  unsigned int value;
  //unsigned int value1;
  //unsigned int value2;
} TNode;
//=============================================================
void RSort_Parallel(TNode *m, unsigned int n)
{
  // Количество задействованных потоков
  unsigned char threads = omp_get_num_procs();
  // Выделяем память под временный массив
  TNode *m_temp = (TNode*)malloc(sizeof(TNode) * n);
  unsigned int *s = (unsigned int*)malloc(sizeof(unsigned int) * 256 * threads);

  #pragma omp parallel num_threads(threads)
  {
    TNode *source = m;
    TNode *dest = m_temp;
    unsigned int l = omp_get_thread_num();
    unsigned int div = n / omp_get_num_threads();
    unsigned int mod = n % omp_get_num_threads();
    unsigned int left_index  = l < mod ? (div + (mod == 0 ? 0 : 1)) * l : n - (omp_get_num_threads() - l) * div;
    unsigned int right_index = left_index + div - (mod > l ? 0 : 1);

    for (unsigned int digit=0; digit< sizeof(m->key); digit++)
    {
      unsigned int s_sum[256] = {0};
      unsigned int s0[256] = {0};
      unsigned char *b1 = (unsigned char*)&source[right_index].key;
      unsigned char *b2 = (unsigned char*)&source[left_index].key;
      while (b1 >= b2)
      {
        s0[*(b1+digit)]++;
        b1 -= sizeof(TNode);
      }
      for (unsigned int i=0; i<256; i++)
      {
        s[i+256*l] = s0[i];
      }

      #pragma omp barrier
      for (unsigned int j=0; j<threads; j++)
      {
        for (unsigned int i=0; i<256; i++)
        {
          s_sum[i] += s[i+256*j];
          if (j<l)
          {
            s0[i] += s[i+256*j];
          }
        }
      }

      for (unsigned int i=1; i<256; i++)
      {
        s_sum[i] += s_sum[i-1];
        s0[i] += s_sum[i-1];
      }
      unsigned char *b = (unsigned char*)&source[right_index].key + digit;
      TNode *v1 = &source[right_index];
      TNode *v2 = &source[left_index];
      while (v1 >= v2)
      {
        dest[--s0[*b]] = *v1--;
        b -= sizeof(TNode);
      }
      #pragma omp barrier
      TNode *temp = source;
      source = dest;
      dest = temp;
    }
  }

  // Если ключ структуры однобайтовый, просто копируем в исходный массив
  if (sizeof(m->key)==1)
  {
    memcpy(m, m_temp, n * sizeof(TNode));
  }

  free(s);
  free(m_temp);
}
//=============================================================
int main()
{
  unsigned int n=10000000;

  TNode *m1 = (TNode*)malloc(sizeof(TNode) * n);
  srand(time(NULL));

  for (unsigned int i=0; i<n; i++)
  {
    m1[i].key = rand()*RAND_MAX+rand();
  }

  RSort_Parallel(m1, n);

  free(m1);
  return 0;
}


Надеюсь было интересно.

С уважением, Ваш покорный слуга, Rebuilder.

P.S.

Сравнение алгоритмов по скорости не привожу сознательно. Результаты сравнений, предложения и критика приветствуются.

P.P.S.