Radix sort - alles auspressen





Grüße an alle Algorithmusliebhaber. Ich möchte Ihnen etwas über meine Forschung zum Sortieren im Allgemeinen erzählen und mich mit der Frage des Sortierens von Radix befassen.







Als Entwickler mit langjähriger Erfahrung sah er sich zunehmend einem merkwürdigen Trend in der Softwareentwicklung gegenüber:



Trotz der Entwicklung der Hardware moderner Computer und Verbesserungen der Algorithmen stieg die Leistung des Codes im Allgemeinen nicht nur nicht an, sondern verschlechterte sich an einigen Stellen erheblich.



Ich glaube, dies liegt an der allgemeinen Idee, schnelles Programmieren mit immer leistungsfähigeren Frameworks und Ultra-High-Level- / Skriptsprachen zu bevorzugen. Sprachen wie Ruby oder Python sind unglaublich entwicklerfreundlich. Viel 'syntaktischer Zucker', ich würde sogar 'Honig' sagen, beschleunigt die Entwicklung manchmal, wenn auch nicht um Größenordnungen, aber zu welchem ​​Preis. Als Benutzer ärgere ich mich über den geringen thermischen Wirkungsgrad des Codes. Ich werde einfach über die Menge des verbrauchten Speichers schweigen, aber die wichtigste Ressource der Menschheit ist die Zeit. Es verschwindet spurlos in endlosen Abstraktionen, wird von Code-Analysatoren gestohlen und von intelligenten Müllsammlern gelöscht. Ich möchte nicht in die Vergangenheit zurückkehren, die Vorteile der modernen Entwicklung aufgeben, "teuren" Code schreiben,Ich schlage nur vor, über die mögliche Beseitigung von Leistungsengpässen nachzudenken, wenn dies bei typischen Aufgaben möglich ist. Dies kann häufig durch die Optimierung von Codeabschnitten mit hoher Last erreicht werden.



Das Sortieren kann als eine der grundlegenden Optimierungsaufgaben herausgestellt werden. Das Thema ist so auf und ab erforscht, dass es den Anschein hat, als sei es schwierig, unterwegs etwas Interessantes zu finden. Wir werden es jedoch versuchen.



Wir werden keine kleinen Arrays sortieren (weniger als eine Million Elemente). Selbst wenn dies äußerst ineffizient ist, ist es ziemlich schwierig, die Nachteile zu spüren, da sie durch die Leistung moderner Geräte ausgeglichen werden. Große Datenmengen (Milliarden von Elementen) sind eine andere Sache, die Ausführungsgeschwindigkeit unterscheidet sich stark von einer kompetenten Auswahl eines Algorithmus.



Alle auf Vergleichen basierenden Algorithmen lösen das Sortierproblem im Allgemeinen nicht besser als O (n * Log n). Bei n nimmt der Wirkungsgrad schnell ab und es ist nicht möglich, diese Situation zu ändern. Dieser Trend kann korrigiert werden, indem vergleichsbasierte Methoden aufgegeben werden. Das vielversprechendste für mich ist der Radix-Sortieralgorithmus. Seine Rechenkomplexität ist O (k * n), wobei k die Anzahl der Durchgänge durch das Array ist. Wenn n groß genug ist, aber k im Gegenteil sehr klein ist, gewinnt dieser Algorithmus O (n * Log n).

In der modernen Prozessorarchitektur wird k fast immer auf die Anzahl der Bytes der sortierten Anzahl reduziert. Zum Beispiel für DWord (int) k = 4, was nicht viel ist. Dies ist eine Art "potenzielle Grube" von Berechnungen, die auf mehrere Faktoren zurückzuführen ist:



  • Prozessorregister werden für 8-Bit-Operatoren auf Hardwareebene geschärft
  • Der im Algorithmus verwendete Zählpuffer passt in eine Zeile des L1-Cache - den Prozessor. (256 * 4-Byte-Zahlen)


Sie können versuchen, die Richtigkeit dieser Aussage selbst zu überprüfen. Gegenwärtig ist jedoch die Division durch Bit die beste Option. Ich schließe nicht aus, dass die Option zum Teilen entlang der 2-Byte-Grenze rentabler wird, wenn der L1-Cache von Prozessoren auf 256 KB wächst.



Eine effiziente Implementierung der Sortierung ist nicht nur ein Algorithmus, sondern auch ein heikles technisches Problem bei der Optimierung von Code.



In dieser Lösung besteht der Algorithmus aus mehreren Stufen:



  1. , ,
  2. ,
  3. (LSD),


Wir wenden den LSD-Algorithmus (zumindest in meiner Version) aufgrund einer reibungsloseren Verarbeitung mit verschiedenen Schwankungen der Eingabedaten schneller an.

Das ursprüngliche vollständig sortierte Array ist der schlechteste Fall für den Algorithmus, da die Daten weiterhin vollständig sortiert werden. Andererseits ist die Radix-Sortierung bei zufälligen oder gemischten Daten äußerst effektiv.



Das Sortieren eines einfachen Array von Zahlen ist selten. Normalerweise benötigen Sie ein Wörterbuch der Form: Schlüsselwert, wobei der Wert ein Index oder ein Zeiger sein kann.

Zur Universalisierung wenden wir eine Struktur der Form an:



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;

      
      





Je kleiner die Bitterkeit des Schlüssels ist, desto schneller arbeitet der Algorithmus natürlich. Da der Algorithmus nicht mit Zeigern auf eine Struktur arbeitet, sondern diese tatsächlich im Speicher steuert. Mit einem Minimum an Feldern erreichen wir eine hohe Geschwindigkeit. Mit zunehmendem Volumen der Datenfelder der Struktur sinkt die Effizienz dramatisch.



Ich habe bereits zuvor eine Notiz über die Implementierung der Radix-Sortierung in Pascal geschrieben, aber die "Trennung" von Programmiersprachen gewinnt unter den Stammgästen dieser Ressource an beispiellosen Raten. Aus diesem Grund habe ich beschlossen, den Code für diesen Artikel in 'si' teilweise als 'korrekter' umzuschreiben. eine gemeinsame Sprache in der Gemeinschaft der Programmierer (obwohl am Ausgang die Generierung von Assembler-Code mit Pascal stellenweise fast identisch ist). Beim Zusammenbau empfehle ich die Verwendung des Schlüssels –O2 oder höher.



Im Gegensatz zur vorherigen Implementierung war es unter Verwendung von Zeigern in der Schleifenadressierung anstelle von bitweisen Verschiebungsoperationen möglich, eine noch größere Erhöhung der Algorithmusgeschwindigkeit um 1-2% zu erzielen.



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.   

      
      







Ich habe sehr lange über diesen Code meditiert, konnte aber nicht besser schreiben. Vielleicht können Sie mir sagen, wie ich diesen Code schneller machen kann.



Und wenn Sie etwas schneller wollen?



Der nächste logische Schritt ist, wie ich es gesehen habe, die Verwendung einer Grafikkarte. Im Internet habe ich viele Gründe dafür gesehen, dass die Radix-Sortierung auf vielen Grafikkartenkernen perfekt parallel ist (fast die beste Sortiermethode ist keine Grafikkarte).



Nachdem ich der Versuchung erlegen war, zusätzliche Leistung zu erzielen, implementierte ich mehrere mir bekannte Sortieroptionen mit OpenCL. Leider habe ich nicht die Möglichkeit, die Implementierung auf Top-End-Grafikkarten zu überprüfen, aber auf meiner GEFORCE GTX 750 TI ging der Algorithmus durch die Single-Threaded-Implementierung auf der CPU verloren, da die Daten an die Grafikkarte gesendet und dann zurückgenommen werden mussten. Wenn Sie die Daten nicht im Bus hin und her laufen lassen, ist die Geschwindigkeit akzeptabel, aber immer noch kein Springbrunnen. Es gibt noch eine Bemerkung. In OpenCL sind die Ausführungsthreads nicht synchron (Arbeitsgruppen werden meines Wissens in einer beliebigen Reihenfolge ausgeführt, dies ist in CUDA nicht der Fall, richtig, wer weiß), was Sie in diesem Fall daran hindert, effizienteren Code zu schreiben.



Code aus der Serie: Ist es möglich, einen Trolleybus zu erstellen ... Verarbeitung in OpenCL bei Delphi ... aber warum?
Ich war zu faul, um in 'C' umzuschreiben, ich poste es so wie es ist.

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.  

      
      







Es gibt noch eine ungetestete Option - Multithreading. Ich habe bereits nur Bare C und die OpenMP-Bibliothek verwendet und mich entschlossen herauszufinden, welche Auswirkungen die Verwendung mehrerer CPU-Kerne haben würde.



Ursprünglich bestand die Idee darin, das ursprüngliche Array in gleiche Teile aufzuteilen, sie in separate Streams zu übertragen und sie dann zusammenzuführen (Sortierung zusammenführen). Die Sortierung verlief gut, aber die Zusammenführung verlangsamte die gesamte Struktur erheblich. Jedes Kleben entspricht einem zusätzlichen Durchgang durch das Array. Der Effekt war schlimmer als das Arbeiten in einem Thread. Die Implementierung wurde als nicht praktikabel abgelehnt.



Infolgedessen habe ich eine parallele Sortierung angewendet, die der auf der GPU verwendeten sehr ähnlich ist. Bei ihr lief alles viel besser.



Merkmale der Parallelverarbeitung:

In der aktuellen Implementierung hat er Probleme mit der Thread-Synchronisation so weit wie möglich vermieden (atomare Funktionen werden auch nicht verwendet, da verschiedene Threads unterschiedliche, wenn auch manchmal benachbarte Speicherblöcke lesen). Threads sind keine freie Sache, der Prozessor verbringt wertvolle Mikrosekunden mit ihrer Erstellung und Synchronisation. Indem Sie die Barrieren auf ein Minimum beschränken, können Sie ein wenig sparen. Da alle Threads denselben L3-Cache-Speicher und RAM verwenden, ist der Gesamtgewinn des Algorithmus aufgrund des Amdahlschen Gesetzes leider nicht so signifikant, da die Anzahl der Threads zunimmt.



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

      
      







Ich hoffe es war interessant.



Mit freundlichen Grüßen Mit freundlichen Grüßen Rebuilder.



PS Ein

Vergleich von Algorithmen hinsichtlich der Geschwindigkeit ist bewusst nicht. Vergleichsergebnisse, Vorschläge und Kritik sind willkommen.



PPS





All Articles