Einführung in Paralleles Programmieren

Startbeispiel

 int feld[32],zielfeld[32];
 for(i=0;i<32;i++)
  {
   feld[i]=i+1; //zu quadrierende Zahlen ins Feld einfuellen
  }
 quadrieren(feld,zielfeld,32); //alle 32 Zahlen quadrieren
 for(i=0;i<32;i++) //Ergebnisse anzeigen
  {
   printf(" %d",zielfeld[i]);
  }

Serielle Version:

void quadrieren(int *feld,int *ziel,int n)
{
 for(int i=0;i<n;i++)
  {
   ziel[i] = feld[i]*feld[i];
  }
}

Version mit Threads auf der CPU:

#include <thread>

void quadr(int *feld,int *ziel,int i)
{
 ziel[i] = feld[i]*feld[i];
}

void quadrieren(int *feld,int *ziel,int n)
{
 std::thread tr[n-1];
 for(int i=0;i<n-1;i++)
  {
   //fuer jeden Thread das Unterprogramm quadr() starten:
   tr[i] = std::thread(quadr,feld,ziel,i);
  }
 quadr(feld,ziel,n-1); //Haupt-Thread startet quadr() direkt
 for(int i=0;i<n-1;i++)
  {
   tr[i].join(); //auf Beenden aller Threads warten
  }
}

Wenn man eine CPU mit mindestens 32 Kernen hat, läuft das Programm somit etwa 32 mal schneller.
Bei weniger Kernen funktioniert das Programm trotzdem, es ist dann aber etwas langsamer, weil das Betriebssystem dann zwischen den verschiedenen Threads wechseln muss.
Man kann das Programm aber auch etwas anpassen, so dass nur soviele Threads verwendet werden wie Prozessorkerne vorhanden sind. Ein entsprechendes Beispiel etwas später.

Version auf GPU mit OpenCL:

__kernel
void quadr(global int *feld,global int *ziel)
{
 int i = get_global_id(0);
 ziel[i] = feld[i]*feld[i];
}

  kernel_quadr.setArg(0,d_feld);
  kernel_quadr.setArg(1,d_ziel);
  queue.enqueueNDRangeKernel(kernel_quadr,cl::NullRange,cl::NDRange(32),cl::NullRange);
  queue.finish(); //auf Beenden aller Threads warten

Version auf GPU mit CUDA:

__global__
void cuda_quadr(int *feld,int *ziel)
{
 int i = threadIdx.x;
 ziel[i] = feld[i]*feld[i];
}

  int nBlocks=1, threadsPerBlock=32;
  cuda_quadr<<<nBlocks,threadsPerBlock>>>(d_feld,d_ziel);

Kopieren der Daten zur Grafikkarte

Leider kann es aber nicht direkt auf den Hauptspeicher (RAM) des Computers zugreifen. Wir müssen deshalb die Daten erst in den Grafikkarten-Speicher kopieren und nach der Berechnung dann die Resultate wieder zurückkopieren.

Im Prinzip muss es nicht unbedingt eine Grafikkarte sein, es wird deshalb meistens von einem Device gesprochen.

Zum Device kopieren geht mit CUDA so:

 int *d_feld; //Device-Variable (Adresse ist Zeiger auf int)
 int size=n*sizeof(int);
 cudaMalloc((void**)&d_feld,size);
 cudaMemcpy(d_feld,feld,size,cudaMemcpyHostToDevice);

 cl::Buffer dfeld; //Device-Variable (Adresse ist ein OpenCL-Objekt)
 int size=n*sizeof(int);
 dfeld = cl::Buffer(context,CL_MEM_READ_WRITE,size);
 queue=cl::CommandQueue(context,default_device);
 queue.enqueueWriteBuffer(dfeld,CL_TRUE,0,size,feld);

Vom Device zurückkopieren geht mit OpenCL so:

 queue.enqueueReadBuffer(dziel,CL_TRUE,0,n*sizeof(int),ziel);

 cudaMemcpy(ziel,d_ziel,n*sizeof(int),cudaMemcpyDeviceToHost);

Übersetzen (compilieren) des Kernels

Die Besonderheit von CUDA ist, dass der Kernel mit einer speziellen Syntax aufgerufen werden muss. Dazu muss man das Developer-Kit von Nvidia installieren und den darin enthaltenen Compiler "nvcc" verwenden.

Das gesamte Beispiel1 habe ich in 4 Teile aufgeteilt:
beispiel1.cc              Hauptprogramm zum beide Varianten aufrufen
beispiel1_cuda.cu     Programmteile fuer CUDA (mit nvcc übersetzen)
beispiel1_opencl.cc   Programmteil fuer OpenCL
beispiel1_kernel.cc    Kernel fuer OpenCL
Dann noch "beispiel1.h" zur Vordeklaration einiger Funktionen und "makefile" damit das Ganze bequem mit "make" compiliert werden kann.

Die interressantesten Teile des Programms sehen somit so aus:

// beispiel1.cc
...
 else if(variante==3)
  {printf("Auf Grafikkarte mit OpenCL gerechnet:\n");
   copy_to_device(feld,32);
   bool ok=opencl_kernel_compilieren("beispiel1_kernel.cc");
   if(!ok) {printf("Fehler3: opencl_kernel_compilieren() misslungen.\n"); exit(1);}
   quadrieren_opencl(feld,zielfeld,32);
   device_speicher_freigeben();
  }
 else //if(variante==4)
  {printf("Auf Grafikkarte mit CUDA gerechnet:\n");
   copy_to_device_cuda(feld,32);
   quadrieren_cuda(feld,zielfeld,32);
   device_speicher_freigeben_cuda();
  }
...

// beispiel1_opencl.cc
...
extern bool copy_to_device(const int *feld,int n)
{
 opencl_init();
 size1 = sizeof(int)*n;
 dfeld = cl::Buffer(context,CL_MEM_READ_WRITE,size1);
 dziel = cl::Buffer(context,CL_MEM_READ_WRITE,size1);
 queue=cl::CommandQueue(context,default_device);
 queue.enqueueWriteBuffer(dfeld,CL_TRUE,0,size1,feld);
 return true;
}

extern bool copy_from_device(int *ziel,int n)
{
 queue.enqueueReadBuffer(dziel,CL_TRUE,0,sizeof(int)*n,ziel);
 return true;
}

extern bool opencl_kernel_compilieren(const char *kernel_name)
{
 if(size1==0)
  {printf("Fehler1: kein Speicher auf dem Device reserviert.\n"); return false;}
 char *quellcode;
 int len=filelaenge(kernel_name);
 if(len==0)
  {printf("Fehler2: \"%s\" nicht gefunden\n",kernel_name); return false;}
 quellcode=new char[len]; if(quellcode==NULL) return false;
 datei_einlesen(kernel_name,quellcode);
 std::string code=quellcode;
 delete[] quellcode;
 cl::Program::Sources sources;
 sources.push_back({code.c_str(),code.length()});
 cl::Program program(context,sources);
 if(program.build({default_device})!=CL_SUCCESS) return false;
 kernel_quadr=cl::Kernel(program,"quadr");
 return true;
}

extern void quadrieren_opencl(const int *feld,int *ziel,int n)
{
 if(size1==0) copy_to_device(feld,n);
 kernel_quadr.setArg(0,dfeld);
 kernel_quadr.setArg(1,dziel);
 queue.enqueueNDRangeKernel(kernel_quadr,cl::NullRange,cl::NDRange(n),cl::NullRange);
 queue.finish(); //auf Beenden aller Threads warten
 copy_from_device(ziel,n);
}

extern void device_speicher_freigeben()
{
 //nicht noetig?
}
...

// beispiel1_cuda.cu
...
int *d_feld=NULL, *d_ziel=NULL;

bool copy_to_device_cuda(const int *feld,int n)
{
 int size=n*sizeof(int);
 cudaMalloc((void**)&d_feld, size);
 cudaMalloc((void**)&d_ziel, size); if(d_ziel==NULL) return false;
 cudaMemcpy(d_feld,feld,size,cudaMemcpyHostToDevice);
 return true;
}

bool copy_from_device_cuda(int *ziel,int n)
{
 cudaMemcpy(ziel,d_ziel,n*sizeof(int),cudaMemcpyDeviceToHost);
 return true;
}

extern void quadrieren_cuda(const int *feld,int *ziel,int n)
{
 if(d_ziel==NULL)
  {
   if(!copy_to_device_cuda(feld,n))
     {printf("Fehler: copy_to_device_cuda() misslungen\n"); return;}
  }
 int nBlocks=1, threadsPerBlock=32;
 cuda_quadr<<<nBlocks,threadsPerBlock>>>(d_feld,d_ziel);
 cudaDeviceSynchronize(); //auf Beenden aller Threads warten
 copy_from_device_cuda(ziel,n);
}

extern void device_speicher_freigeben_cuda()
{
 cudaFree(d_ziel); d_ziel=NULL;
 cudaFree(d_feld); d_feld=NULL;
}
...

Compilieren des gesamten Programms

Unter Linux (Ubuntu) braucht man folgende Pakete:
build-essential
opencl-headers
opencl-c-headers ?
opencl-cpp-headers ?
ocl-icd-*
ocl-icd-opencl-dev
ocl-icd-libopencl1
beignet (für Intel-Grafikkarten)
beignet-devel (für Intel-Grafikkarten)
clinfo
nvidia-opencl-dev
nvidia-cuda-dev
nvidia-cuda-doc
nvidia-cuda-gdb
nvidia-cuda-toolkit
libcudart10.1
nvidia-390 ?
nvidia-utils-390 ?

Versionsnummern können ändern. Bei älteren Grafikkarten sind eventuell auch ältere Pakete erforderlich.
Das "CUDA Toolkit" und der zur Grafikkarte passende Nvidia-Treiber muss man direkt von der Nvidia-Homepage runterladen.
Eventuell werden dabei einige der obigen Pakete schon installiert?

Wenn alle Pakete installiert sind, kann das Programm so compiliert und laufen gelassen werden:

 tar zxvf beispiel1.tar.gz
 cd beispiel1/
 make
 ./beispiel1 2
 ./beispiel1 3
 ./beispiel1 4

Geschwindigkeits-Tests

#define N 1000000
 int *feld,*zielfeld;
 feld = new int[N];
 zielfeld = new int[N];
 for(i=0;i<N;i++)
  {
   feld[i] = random()&0xFFFF;
   zielfeld[i]=0;
  }

 delete[] feld;
 delete[] zielfeld;

Der entsprechende Programmteil sieht dann so aus:

void quadr(int *feld,int *ziel,int i0,int nt,int n)
{
 for(int i=i0;i<n;i+=nt)
   ziel[i] = feld[i]*feld[i];
}

void quadrieren2(int *feld,int *ziel,int n)
{
 //Anzahl verfuegbare Threads:
 int nt = std::thread::hardware_concurrency();
 std::thread tr[nt-1];
 for(int i=0;i<nt-1;i++)
   tr[i]=std::thread(quadr,feld,ziel,i,nt,n);
 quadr(feld,ziel,nt-1,nt,n);
 for(int i=0;i<nt-1;i++)
   tr[i].join();
}

Die Variante mit OpenCL können wir so lassen.
Wenn der Compiler genügend intelligent ist, sorgt er dafür dass die Arbeit auf die verfügbare Anzahl Threads aufgeteilt wird.
Aus Programmierersicht sieht es so aus, wie wenn 1000000 Threads parallel gerechnet würden.

Die Variante mit CUDA müssen wir noch etwas anpassen.
Ein Block kann nur eine bestimmte Anzahl Threads wirklich parallel ausführen. Ich glaube auf den bisher getesteten Grafikkarten sind das 4 bis 32. Man kann zwar bis zu 1024 angeben, aber die werden dann nicht wirklich parallel gerechnet.
Allerdings sind eine grosse Anzahl Blöcke pro Grid vorhanden. Bis zu 65535, wobei soviele auch nicht wirklich physikalisch vorhanden sind, womit dann einige Blöcke doch wieder nacheinander abgearbeitet werden.
Die Anzahl Blöcke die wir brauchen, so dass theoretisch alles parallel gerechnet wird, können wir leicht ausrechnen.
Wenn wir Anzahl Threads pro Block auf 32 setzen, dann sind das: N/32 auf ganze Zahl aufgerundet.
Eventuell bekommen dann einige (wenige) Threads eine id, die zu gross für unser Zahlenfeld ist. Das muss dann im Kernel abgefangen werden.
Die entsprechenden Programmteile in beispiel1_cuda.cu:

__global__
void cuda_quadr(const int *feld,int *ziel,int n)
{
 int i = blockDim.x*blockIdx.x + threadIdx.x;
 if(i<n)
  {
   ziel[i] = feld[i]*feld[i];
  }
}
...
extern void quadrieren_cuda(const int *feld,int *ziel,int n)
{
 if(d_ziel==NULL)
  {
   if(!copy_to_device_cuda(feld,n))
     {printf("Fehler1: copy_to_device_cuda() misslungen\n"); return;}
  }
 int threadsPerBlock=32; //bis 1024 erlaubt
 int nBlocks=divup(n,threadsPerBlock); //bis 65535 erlaubt
 //(divup(a,b) ist aufgerundete Division a/b)
 cuda_quadr<<<nBlocks,threadsPerBlock>>>(d_feld,d_ziel,n);
 cudaDeviceSynchronize(); //auf Beenden aller Threads warten
 errcheck("Fehler2: Kernelaufruf misslungen");
 copy_from_device_cuda(ziel,n);
}

Wenn wir auf mehr als 65535 Anzahl Blöcke kommen, müssen wir entweder die Anzahl Threads pro Block erhöhen, oder wir machen es ähnlich wie bei den CPU-Threads.
Die gesamte Anzahl der Threads können wir dabei wieder durch fest definierte Register erhalten.
Der Kernel sieht dann so aus:

__global__
void cuda_quadr(const int *feld,int *ziel,int n)
{
 int nt = blockDim.x*gridDim.x; //Anzahl Threads
 int i0 = blockDim.x*blockIdx.x + threadIdx.x;
 for(int i=i0; i<n; i+=nt)
  {
   ziel[i] = feld[i]*feld[i];
  }
}

 int threadsPerBlock=32;
 int nBlocks=1024; //bis 65535 erlaubt
 cuda_quadr<<<nBlocks,threadsPerBlock>>>(d_feld,d_ziel,n);

Für Zeitmessungen habe ich noch eine Option -q ins Hauptprogramm eingebaut um die Überprüfung der Ergebnisse wegzulassen.

Das so erweiterte erste Beispiel kann hier heruntergeladen werden:
beispiel1b.tar.gz
Übersetzen wie gehabt.
Zum laufen lassen mit Zeitmessung so aufrufen:
time ./beispiel1 4
Und mit Zeitmessung ohne Resultate zu prüfen so aufrufen:
time ./beispiel1 -q 4

Zu schnell für Zeitmessung!?

Idee1:

Idee2:

Idee3:

Alle diese Ideen waren nicht ideal. Ich stelle dann im Teil2 eine bessere Methode vor.

Fortsetzungen:

Anmerkungen:

Grid: