Radix Sort - tout extraire





Salutations à tous les amateurs d'algorithmes. Je veux vous parler de mes recherches sur le tri en général et approfondir la réflexion sur le tri à la base.







En tant que développeur avec de nombreuses années d'expérience, il a de plus en plus commencé à faire face à une étrange tendance dans le développement de logiciels:



malgré le développement du matériel des ordinateurs modernes et les améliorations des algorithmes, en général, les performances du code n'ont pas seulement augmenté, mais dans certains endroits, elles se sont à peu prÚs dégradées.



Je pense que cela est dĂ» Ă  l'idĂ©e gĂ©nĂ©rale de donner la prĂ©fĂ©rence Ă  la programmation rapide utilisant des frameworks de plus en plus puissants et des langages de script / de trĂšs haut niveau. Des langages comme Ruby ou Python sont incroyablement conviviaux pour les dĂ©veloppeurs. Beaucoup de «sucre syntaxique», je dirais mĂȘme «miel», accĂ©lĂšrent parfois le dĂ©veloppement, sinon des ordres de grandeur, mais Ă  quel prix. En tant qu'utilisateur, je suis ennuyĂ© par la faible efficacitĂ© thermique du code, je vais simplement garder le silence sur la quantitĂ© de mĂ©moire consommĂ©e, mais la principale ressource de l'humanitĂ© est le temps. Il disparaĂźt sans laisser de trace dans des abstractions sans fin, est volĂ© par des analyseurs de code, effacĂ© par des garbage collors intelligents. Je n’ai pas envie de revenir dans le passĂ©, d’abandonner les avantages du dĂ©veloppement moderne, d’écrire du code «coĂ»teux»,Je propose simplement de rĂ©flĂ©chir Ă  l'Ă©limination possible des goulots d'Ă©tranglement de performance lorsque cela est possible dans les tĂąches typiques. Cela peut souvent ĂȘtre rĂ©alisĂ© en optimisant les sections de code Ă  forte charge.



Le tri peut ĂȘtre distinguĂ© comme l'une des tĂąches d'optimisation de base. Le sujet est tellement explorĂ©, de haut en bas, qu'il semblerait qu'il soit difficile de trouver quelque chose d'intĂ©ressant en cours de route. Cependant, nous essaierons.



Nous ne trierons pas les petits tableaux (moins d'un million d'Ă©lĂ©ments). MĂȘme s'il est extrĂȘmement inefficace de le faire, il est assez difficile de ressentir les baisses, car elles sont nivelĂ©es par les performances des Ă©quipements modernes. De grandes quantitĂ©s de donnĂ©es (des milliards d'Ă©lĂ©ments) sont une autre question; la vitesse d'exĂ©cution varie considĂ©rablement d'une sĂ©lection compĂ©tente d'un algorithme.



Tous les algorithmes basĂ©s sur des comparaisons rĂ©solvent gĂ©nĂ©ralement le problĂšme de tri pas mieux que O (n * Log n). Au grand n, l'efficacitĂ© diminue rapidement et il n'est pas possible de changer cette situation. Cette tendance peut ĂȘtre corrigĂ©e en abandonnant les mĂ©thodes basĂ©es sur la comparaison. Le plus prometteur pour moi est l'algorithme de tri Radix. Sa complexitĂ© de calcul est O (k * n), oĂč k est le nombre de passages dans le tableau. Si n est assez grand, mais k, au contraire, est trĂšs petit, alors cet algorithme l'emporte sur O (n * Log n).

Dans l'architecture de processeur moderne, k est presque toujours réduit au nombre d'octets du nombre trié. Par exemple, pour DWord (int) k = 4, ce qui n'est pas beaucoup. Il s'agit d'une sorte de «gouffre potentiel» de calculs, qui est dû à plusieurs facteurs:



  • Les registres du processeur sont affinĂ©s pour les opĂ©rateurs 8 bits au niveau matĂ©riel
  • Le tampon de comptage utilisĂ© dans l'algorithme s'inscrit dans une ligne du cache L1 - le processeur. (256 * nombres de 4 octets)


Vous pouvez essayer de vĂ©rifier vous-mĂȘme la vĂ©racitĂ© de cette dĂ©claration. Cependant, Ă  l'heure actuelle, la division par bit est la meilleure option. Je n'exclus pas que lorsque le cache L1 des processeurs passe Ă  256 Ko, l'option de diviser le long de la limite de 2 octets deviendra plus rentable.



Une mise en Ɠuvre efficace du tri n'est pas seulement un algorithme, mais aussi un problĂšme d'ingĂ©nierie dĂ©licat pour optimiser le code.



Dans cette solution, l'algorithme se compose de plusieurs Ă©tapes:



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


Nous appliquons l'algorithme LSD plus rapidement (du moins dans ma version) en raison d'un traitement plus fluide avec diverses fluctuations des données d'entrée.

Le tableau d'origine entiĂšrement triĂ© est le pire des cas pour l'algorithme, car les donnĂ©es seront toujours entiĂšrement triĂ©es. En revanche, le tri Radix est extrĂȘmement efficace sur des donnĂ©es alĂ©atoires ou mixtes.



Le tri d'un simple tableau de nombres est rare, gĂ©nĂ©ralement vous avez besoin d'un dictionnaire de la forme: clĂ© - valeur, oĂč la valeur peut ĂȘtre un index ou un pointeur.

Pour l'universalisation, nous appliquerons une structure de la forme:



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;

      
      





Naturellement, plus le bitness de la clé est petit, plus l'algorithme fonctionne rapidement. Puisque l'algorithme ne fonctionne pas avec des pointeurs vers une structure, mais le pilote en fait en mémoire. Avec un minimum de champs, nous obtenons une vitesse élevée. Avec une augmentation du volume des champs de données de la structure, l'efficacité diminue considérablement.



Auparavant, j'avais dĂ©jĂ  Ă©crit une note avec l'implĂ©mentation du tri Radix en Pascal, cependant, la «sĂ©grĂ©gation» des langages de programmation gagne des taux sans prĂ©cĂ©dent chez les habituĂ©s de cette ressource. Par consĂ©quent, j'ai dĂ©cidĂ© de rĂ©Ă©crire partiellement le code de cet article dans 'si' comme plus 'correct'un langage commun dans la communautĂ© des programmeurs (bien que la sortie de la gĂ©nĂ©ration de code assembleur avec Pascal soit presque identique par endroits). Lors de l'assemblage, je recommande d'utiliser la clĂ© –O2 ou supĂ©rieure.



Contrairement à l'implémentation précédente, en utilisant des pointeurs dans l'adressage en boucle au lieu d'opérations de décalage au niveau du bit, nous avons réussi à obtenir une augmentation encore plus grande de 1 à 2% de la vitesse de l'algorithme.



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.   

      
      







J'ai mĂ©ditĂ© sur ce code pendant trĂšs longtemps, mais je ne pouvais pas mieux Ă©crire. Peut-ĂȘtre pouvez-vous me dire comment rendre ce code plus rapide.



Et si vous voulez un peu plus vite?



La prochaine Ă©tape logique, comme je l'ai vu, est d'utiliser une carte vidĂ©o. Sur Internet, j'ai vu beaucoup de raisonnements selon lesquels le tri Radix est parfaitement parallĂšle sur de nombreux cƓurs de carte vidĂ©o (presque la meilleure mĂ©thode de tri n'est pas une carte vidĂ©o).



Ayant succombĂ© Ă  la tentation d'obtenir des performances supplĂ©mentaires, j'ai implĂ©mentĂ© plusieurs options de tri en utilisant OpenCL que je connais. Malheureusement, je n'ai pas l'occasion de vĂ©rifier l'implĂ©mentation sur les cartes vidĂ©o haut de gamme, mais sur ma GEFORCE GTX 750 TI, l'algorithme a perdu Ă  l'implĂ©mentation monothread sur le CPU, du fait que les donnĂ©es devaient ĂȘtre envoyĂ©es Ă  la carte vidĂ©o puis reprises. Si vous n’exĂ©cutiez pas les donnĂ©es dans les deux sens dans le bus, la vitesse serait acceptable, mais pas encore une fontaine. Il y a encore une remarque. Dans OpenCL, les threads d'exĂ©cution ne sont pas synchrones (les groupes de travail sont exĂ©cutĂ©s dans un ordre arbitraire, pour autant que je sache, ce n'est pas le cas dans CUDA, correct, qui sait), ce qui vous empĂȘche d'Ă©crire du code plus efficace dans ce cas.



Code de la série: est-il possible de créer un trolleybus ... traitement en OpenCL chez Delphi ... mais pourquoi?
J'étais trop paresseux pour réécrire en «C», je le poste tel quel.

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.  

      
      







Il reste une autre option non testĂ©e: le multithreading. N'utilisant dĂ©jĂ  que Bare C et la bibliothĂšque OpenMP, j'ai dĂ©cidĂ© de savoir quel effet aurait l'utilisation de plusieurs cƓurs de processeur.



Au dĂ©part, l'idĂ©e Ă©tait de diviser le tableau d'origine en parties Ă©gales, de les transfĂ©rer dans des flux sĂ©parĂ©s, puis de les fusionner (tri par fusion). Le tri s'est bien passĂ©, mais la fusion a beaucoup ralenti l'ensemble de la structure, chaque collage Ă©quivaut Ă  un passage supplĂ©mentaire Ă  travers le tableau. L'effet Ă©tait pire que de travailler dans un fil. RejetĂ© la mise en Ɠuvre car pas pratique.



En conséquence, j'ai appliqué un tri parallÚle trÚs similaire à celui utilisé sur le GPU. Avec elle, tout s'est bien mieux passé.



Caractéristiques du traitement parallÚle:

Dans l'implĂ©mentation actuelle, il a Ă©vitĂ© autant que possible les problĂšmes de synchronisation des threads (les fonctions atomiques ne sont pas non plus utilisĂ©es en raison du fait que diffĂ©rents threads lisent des blocs de mĂ©moire diffĂ©rents, bien que parfois voisins). Les threads ne sont pas gratuits, le processeur consacre de prĂ©cieuses microsecondes Ă  leur crĂ©ation et Ă  leur synchronisation. En minimisant les barriĂšres, vous pouvez Ă©conomiser un peu. Malheureusement, comme tous les threads utilisent la mĂȘme mĂ©moire cache L3 et la mĂȘme RAM, le gain global de l'algorithme n'est pas si important en raison de la loi d' Amdahl , avec une augmentation du nombre de threads.



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

      
      







J'espÚre que c'était intéressant.



Meilleures salutations, Bien Ă  vous, Rebuilder.



PS La

comparaison des algorithmes en termes de vitesse ne l'est volontairement pas. Les résultats de comparaison, les suggestions et les critiques sont les bienvenus.



PPS





All Articles