✉️ 🈴 🤶🏿 Radixsort-すべてを絞り出します 🧑🏼‍🤝‍🧑🏻 🐵 🖐🏾

すべてのアルゴリズム愛好家へのご挨拶。ソーティング全般に関する私の研究についてお話しし、ラディックスのソーティングの検討について掘り下げたいと思います。

長年の経験を持つ開発者として、彼はますますソフトウェア開発の奇妙な傾向に直面し始めました。

現代のコンピューターのハードウェアの開発とアルゴリズムの改善にもかかわらず、一般に、コードのパフォーマンスは向上しなかっただけでなく、場所によってはかなり低下しました。

これは、ますます強力なフレームワークと超高レベル/スクリプト言語を使用した高速プログラミングを優先するという一般的な考え方によるものだと思います。 RubyやPythonなどの言語は、開発者にとって非常に使いやすい言語です。多くの「統語的砂糖」、私は「ハニー」とさえ言うでしょう、時には、規模ではないにしても、いくらのコストで開発をスピードアップします。ユーザーとして、私はコードの熱効率が低いことに悩まされています。消費されるメモリの量については黙っていますが、人類の主なリソースは時間です。それは無限の抽象化で痕跡なしに消え、コードアナライザーによって盗まれ、スマートなガベージコレクターによってクリアされます。私は過去に戻り、現代の開発の恩恵を放棄し、「高価な」コードを書き、私は、典型的なタスクで可能な場合、パフォーマンスのボトルネックを排除できる可能性について考えることを提案します。これは多くの場合、コードの高負荷セクションを最適化することで実現できます。

並べ替えは、基本的な最適化タスクの1つとして選択できます。トピックは上下に非常に探索されているため、途中で何か面白いものを見つけるのは難しいように思われます。ただし、試してみます。

小さな配列（100万未満の要素）はソートしません。これを行うのは非常に非効率的ですが、ドローダウンは最新の機器のパフォーマンスによって平準化されているため、ドローダウンを感じることは非常に困難です。大量のデータ（数十億の要素）は別の問題です。実行速度は、アルゴリズムの適切な選択によって大きく異なります。

比較に基づくすべてのアルゴリズムは、一般に、O（n * Log n）よりも優れたソート問題を解決しません。 nが大きいと、効率が急激に低下し、この状況を変えることはできません。この傾向は、比較ベースの方法を放棄することで修正できます。私にとって最も有望なのは、Radixソートアルゴリズムです。その計算の複雑さはO（k * n）です。ここで、kは配列を通過するパスの数です。 nが十分に大きいが、逆にkが非常に小さい場合、このアルゴリズムはO（n * Log n）に勝ちます。

最新のプロセッサアーキテクチャでは、kはほとんどの場合、ソートされた数のバイト数に削減されます。たとえば、DWord（int）k = 4の場合、これはそれほど多くありません。これは、計算のある種の「潜在的な落とし穴」であり、いくつかの要因によるものです。

プロセッサレジスタは、ハードウェアレベルで8ビットオペレータ向けにシャープ化されています
アルゴリズムで使用されるカウントバッファは、L1キャッシュの1行（プロセッサ）に収まります。（256 * 4バイトの数値）

このステートメントの真実を自分で確認してみてください。ただし、現時点では、ビットで割るのが最善のオプションです。プロセッサのL1キャッシュが256KBに増加すると、2バイトの境界に沿って分割するオプションの収益性が高まることを排除しません。

ソートの効率的な実装は、アルゴリズムであるだけでなく、コードを最適化するための微妙なエンジニアリングの問題でもあります。

このソリューションでは、アルゴリズムはいくつかの段階で構成されています。

, ,
,
(LSD),

入力データのさまざまな変動による処理がスムーズになるため、LSDアルゴリズムを（少なくとも私のバージョンでは）より高速なものとして使用します。

データは引き続き完全にソートされるため、元の完全にソートされた配列はアルゴリズムの最悪のケースです。一方、Radixソートは、ランダムデータまたは混合データに対して非常に効果的です。

単純な数値の配列をソートすることはまれです。通常、次の形式のディクショナリが必要です。key-value。値はインデックスまたはポインタです。

普遍化のために、私たちはフォームの構造を使用します

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;

当然、キーのビット数が小さいほど、アルゴリズムの動作は速くなります。アルゴリズムは構造へのポインタでは機能しませんが、実際にはメモリ内で駆動するためです。最小限のフィールドで、高速になります。構造のデータフィールドの量が増えると、効率が劇的に低下します。

以前、PascalでのRadixソートの実装に関するメモをすでに書きましたが、プログラミング言語の「分離」は、このリソースの常連の間で前例のない速度を獲得しています。したがって、この記事のコードを「si」で「より~~正しい」もの~~として部分的に書き直すことにしました。プログラマーのコミュニティで共通の言語（ただし、Pascalを使用したアセンブラーコード生成の出力は場所によってほぼ同じです）。組み立てるときは、–O2以上のキーを使用することをお勧めします。

以前の実装とは対照的に、ビット単位のシフト操作の代わりにループアドレッシングでポインターを使用すると、アルゴリズム速度をさらに1〜2％向上させることができました。

#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;
}

パスカル

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の並べ替えが多くのビデオカードコアで完全に並列であるという多くの理由を目にしました（ほとんどの場合、最適な並べ替え方法はビデオカードではありません）。

追加のパフォーマンスを得たいという誘惑に負けて、私は知っているOpenCLを使用していくつかのソートオプションを実装しました。残念ながら、トップエンドのビデオカードの実装を確認する機会はありませんが、GEFORCE GTX 750 TIでは、データをビデオカードに送信してから元に戻す必要があるため、CPUのシングルスレッド実装でアルゴリズムが失われました。バス上でデータを前後に実行しなかった場合、速度は許容範囲内ですが、それでも噴水ではありません。もう1つコメントがあります。OpenCLでは、実行スレッドは同期していません（私が知る限り、ワークグループは任意の順序で実行されますが、CUDAの場合はそうではありません）。この場合、より効率的なコードを作成できません。

シリーズのコード：DelphiのOpenCLでトロリーバス処理を作成することは可能ですか...しかし、なぜですか？

怠惰すぎて「C」で書き直すことができなかったので、そのまま投稿します。

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.

テストされていないオプションがもう1つ残っています-マルチスレッドです。すでにベアCとOpenMPライブラリのみを使用しているので、複数のCPUコアを使用するとどのような影響があるかを調べることにしました。

当初のアイデアは、元の配列を等しい部分に分割し、それらを別々のストリームに転送してから、それらをマージすることでした（マージソート）。並べ替えは悪くありませんでしたが、マージによって構造全体の速度が大幅に低下しました。各接着は、アレイを通過する追加のパスに相当します。効果は、1つのスレッドで作業するよりも悪かった。実用的ではないとして実装を拒否しました。

その結果、GPUで使用されているものと非常によく似た並列ソートを適用しました。彼女と一緒に、すべてがはるかに良くなりました。

並列処理の機能：

現在の実装では、スレッド同期の問題を可能な限り回避しました（異なるスレッドが異なるメモリブロックを読み取るため、アトミック関数も使用されません）。スレッドは無料のものではありません。プロセッサは、スレッドの作成と同期に貴重なマイクロ秒を費やします。障壁を最小限に抑えることで、少し節約できます。残念ながら、すべてのスレッドが同じL3キャッシュメモリとRAMを使用するため、スレッド数が増えると、Amdahlの法則により、アルゴリズムの全体的なゲインはそれほど重要ではなくなります。

#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;
}

面白かったと思います。

よろしくお願いします、あなたは本当に、リビルダー。

PS

速度の観点からのアルゴリズムの比較は意図的にではありません。比較結果、提案、批判は大歓迎です。

PPS

Radixsort-すべてを絞り出します

More articles: