Einführung in Paralleles Programmieren 3.Teil

Optimierbares Beispiel

Ich möchte hier aber ein einfacheres Beispiel verwenden.
Es sollen zwei Zahlenfelder miteinander multipliziert werden. Also jeder Eintrag im Feld A soll mit jedem Eintrag im Feld B multipliziert werden. Die Summe (ein Eintrag in A mit allen Einträgen in B multipliziert) soll dann in einem Feld C (mit selber Grösse wie A) gespeichert werden.

Zuerst die serielle Variante:

void hello_calc(const float *A,const float *B,float *C,int X,int Y)
{
 for(int x=0;x<X;x++)
  {
   float sum=0;
   for(int y=0;y<Y;y++)
     sum += A[x]*B[y];
   C[x] = sum;
  }
}

Und jetzt der OpenCL-Kernel:

__kernel void hello(__global const float *A,__global const float *B,
                    __global float *C,const int X,const int Y)
{
 int x=get_global_id(0);
 float sum=0;
 if(x<X)
  {
   for(int y=0;y<Y;y++)
     sum += A[x]*B[y];
   C[x] = sum;
  }
}

Die CUDA-Version wäre praktisch gleich. Ich konzentriere mich in diesem Teil aber auf die OpenCL-Version.
Beim Aufruf können wir uns auf die eindimensionale Variante beschränken:

 const int TPB=1024; //ThreadsPerBlock
 queue.enqueueNDRangeKernel(kernel_hello,cl::NullRange,cl::NDRange(rup(X,TPB)),cl::NDRange(TPB));

Im Hauptprogramm wieder die Stoppuhr-Funktionen von Teil2 eingebaut.
Vergleich von CPU- mit OpenCL-Version gibt je nach verwendeter Grafikkarte einen um etwa Faktor 8 bis 600 besseren Wert. (Faktor 8 auf der CPU würde man mit Threads und Optimierung auch noch erreichen.)

erster Optimierungsversuch

Wir können jetzt also die Daten vom Feld B in den lokalen Speicher laden. Danach werden diese für jeden Eintrag in A benutzt um zu multiplizieren und Summe zu bilden.
Der entsprechende Kernel sieht dann so aus:

#define YMAX 12000

__kernel void hello(__global const float *A,__global const float *B,
                    __global float *C,const int X,const int Y)
{
 int x0=get_local_id(0);
 int x=get_global_id(0);

 __local float Blocal[YMAX]; //lokaler Speicher pro Block (max.48KB)
 int ymax=(YMAX<Y)?YMAX:Y;
 for(int y=x0;y<ymax;y+=1024)
  {
   Blocal[y]=B[y]; //gemeinsam gesamtes B in lokalen Speicher kopieren
  }
 barrier(CLK_LOCAL_MEM_FENCE); //Warten bis alle Threads fertig mit
                               //speichern in lokalen Speicher.
 if(x<X)
  {
   float sum=0;
   for(int y=0;y<ymax;y++)
     sum += A[x]*Blocal[y]; //B nur noch aus lokalem Speicher lesen
   C[x] = sum;
  }
}

In der Schlaufe beim kopieren in den lokalen Speicher arbeiten jetzt die 1024 Threads zusammen. In jedem Schlaufendurchgang werden also 1024 Werte kopiert.
(Der Wert 1024 sollte eigentlich noch durch get_local_size(0) ersetzt werden.)

Die Funktion barrier() wartet bis alle Threads im aktuellen Block soweit sind. Dabei ist es wichtig dass dieser Punkt von allen Threads erreicht wird. Also z.B. zuerst ein if(x<X) zu machen wäre falsch, und könnte ein unerwartetes Verhalten verursachen (vermutlich Absturz oder Programm hängt).

Um die Beschränkung auf 12000 zu umgehen, kann man auch mehrmals jeweils nur einen Teil vom Feld B in den lokalen Speicher laden. Man muss nur aufpassen, dass die Barrieren innerhalb der Schlaufe auch wirklich von jedem Thread erreicht werden.
Wir erweitern also die y-Schlaufe. Im äusseren Teil erhöhen wir jeweils um TPB, lesen dann den entsprechenden Teil vom Feld B ein, und benutzen in der inneren Schlaufe den lokalen Speicher.
Der Kernel sieht dann etwa so aus:

__kernel void hello(__global const float *A,__global const float *B,
		    __global float *C,const int X,const int Y)
{
 int x0=get_local_id(0);
 int x=get_global_id(0);
 const int TPB=get_local_size(0);
 __local float Blocal[1024]; //4KB lokaler Speicher
 float sum=0;
 int y;
 for(y=0;y<Y-TPB;y+=TPB)
  {
   Blocal[x0]=B[y+x0]; //aktuellen Teil von B in lokalen Speicher kopieren
   barrier(CLK_LOCAL_MEM_FENCE); //Warten bis alle Threads fertig
   for(int yb=0;yb<TPB;yb++)
      {
       sum += A[x]*Blocal[yb]; //B nur noch aus lokalem Speicher lesen
      }
   barrier(CLK_LOCAL_MEM_FENCE);
  }
 Blocal[x0]=B[y+x0]; //letzter Teil von B in lokalen Speicher kopieren
 barrier(CLK_LOCAL_MEM_FENCE);
 for(int yb=0;y+yb<Y;yb++)
  {
   sum += A[x]*Blocal[yb]; //B nur noch aus lokalem Speicher lesen
  }
 if(x<X) C[x] = sum;
}

Auf zwei Computern getestet hat diese Optimierung auf dem einen ein Faktor 2 (GT 630 OEM) gebracht, auf dem andern (RTX 2080) praktisch keine Veränderung. (ist auf dieser neueren Grafikkarte der globale Speicher auch schon so schnell?)

2. Optimierungsversuch

In OpenCL ist der Vektortyp float4 in etwa so definiert:

struct float4 {float x,y,z,w;}

Wir benutzen in unserem Kernel also mal float4 für die Felder A und C:

__kernel void hello(__global const float4 *A,__global const float *B,
                    __global float4 *C,const int X,const int Y)
{
 int x=get_global_id(0);
 float4 sum={0,0,0,0};
 if(x<X/4)
  {
   for(int y=0;y<Y;y++)
    {
     sum += B[y]*A[x];
    }
   C[x] = sum;
  }
}

Beim Aufruf sollte auch noch durch 4 dividiert werden (mit Aufrundung):

 queue.enqueueNDRangeKernel(kernel_hello,cl::NullRange,cl::NDRange(rup((X+3)/4,TPB)),cl::NDRange(TPB));

3. Optimierungsversuch

Wenn jetzt aber als Anzahl Threads pro Block statt 1024 nur 64 gesetzt wird, haben wir wieder eine Verbesserung. Ich habe keine richtige Erklärung dafür. Möglicherweise könnte es damit zusammenhängen, dass wenn weniger Threads pro Block verwendet werden, dann mehr Register pro Thread benutzbar sind.

Jedenfalls habe ich als optimale Zahl durch probieren folgende gefunden:
RTX-2080:   64
GT-630-OEM: 128

Bei mir gibt es damit folgende Verbesserungen: 5160 GFLOPS beim RTX-2080, 166 GFLOPS beim GT-630-OEM.
Damit liegt der Wert für den GT-630 leicht über dem Vergleichswert der Matrixmultiplikation mit CUBLAS-Library.

4. Optimierungsversuch

__kernel void hello(__global const float8 *A,__global const float8 *B,
		    __global float8 *C,const int X,const int Y)
{
 int x=get_global_id(0);
 float8 sum={0,0,0,0,0,0,0,0};
 //float8 sum=0;
 if(x<(X+7)/8)
  {
   for(int y=0;y<Y/8;y++)
    {
     sum += B[y].s0*A[x];
     sum += B[y].s1*A[x];
     sum += B[y].s2*A[x];
     sum += B[y].s3*A[x];
     sum += B[y].s4*A[x];
     sum += B[y].s5*A[x];
     sum += B[y].s6*A[x];
     sum += B[y].s7*A[x];
    }
   C[x] = sum;
  }
}

   int Y8=Y/8, Yrest=Y%8;
   if(Yrest!=0)
    {
     sum += B[Y8].s0*A[x];
     if(Yrest>=2) sum += B[Y8].s1*A[x];
     if(Yrest>=3) sum += B[Y8].s2*A[x];
     if(Yrest>=4) sum += B[Y8].s3*A[x];
     if(Yrest>=5) sum += B[Y8].s4*A[x];
     if(Yrest>=6) sum += B[Y8].s5*A[x];
     if(Yrest==7) sum += B[Y8].s6*A[x];
    }
   C[x] = sum;

Wieder die Anzahl Threads pro Block variiert:
RTX-2080:   128
GT-630-OEM: 128...512

Bei mir gibt es damit folgende Verbesserungen: 5594 GFLOPS beim RTX-2080, 214 GFLOPS beim GT-630-OEM.
Also schon sehr nahe am CUBLAS-Vergleichswert, oder sogar darüber.

5. Optimierungsversuch

Hier der entsprechende Kernel:

__kernel void hello(__global const float8 *A,__global const float8 *B,
		    __global float8 *C,const int X,const int Y)
{
 int x0=get_local_id(0);
 const int TPB=get_local_size(0); //Threads Per Block
 int x=get_group_id(0)*TPB+x0; //gleiches Resultat wie get_global_id(0)
 __local float8 Blocal[1024]; //lokaler Speicher (max.48KB)
 float8 sum={0,0,0,0,0,0,0,0};
 const int X8=(X+7)/8;
 const int Y8=Y/8;
 int y;
 for(y=0;y<Y8-TPB;y+=TPB)
  {
   Blocal[x0]=B[y+x0];
   barrier(CLK_LOCAL_MEM_FENCE);
   if(x<X8)
    for(int yb=0;yb<TPB;yb++)
    {
     sum += Blocal[yb].s0*A[x];
     sum += Blocal[yb].s1*A[x];
     sum += Blocal[yb].s2*A[x];
     sum += Blocal[yb].s3*A[x];
     sum += Blocal[yb].s4*A[x];
     sum += Blocal[yb].s5*A[x];
     sum += Blocal[yb].s6*A[x];
     sum += Blocal[yb].s7*A[x];
    }
   barrier(CLK_LOCAL_MEM_FENCE);
  }
 if(y+x0<(Y+7)/8) Blocal[x0]=B[y+x0];
 barrier(CLK_LOCAL_MEM_FENCE);
 if(x<X8)
 {
  for(int yb=0;y+yb<Y8;yb++)
   {
    sum += Blocal[yb].s0*A[x];
    sum += Blocal[yb].s1*A[x];
    sum += Blocal[yb].s2*A[x];
    sum += Blocal[yb].s3*A[x];
    sum += Blocal[yb].s4*A[x];
    sum += Blocal[yb].s5*A[x];
    sum += Blocal[yb].s6*A[x];
    sum += Blocal[yb].s7*A[x];
   }
  switch(Y%8)
   {
    case 7: sum += B[Y8].s6*A[x];
    case 6: sum += B[Y8].s5*A[x];
    case 5: sum += B[Y8].s4*A[x];
    case 4: sum += B[Y8].s3*A[x];
    case 3: sum += B[Y8].s2*A[x];
    case 2: sum += B[Y8].s1*A[x];
    case 1: sum += B[Y8].s0*A[x];
   }
  C[x] = sum;
 }
}

zufällige Optimierungen

Beispiel3 download

Es ist so voreingestellt, dass es beim Compilieren Informationen (für Fehlersuche) anzeigen sollte und nebenbei noch ein tmp.ptx (Assemblerähnlicher Code) speichert.
Mehrmals den gleichen Kernel compiliert gibt die Informationen aber nur beim ersten mal (Bug im OpenCL?).
Das Hauptprogramm heisst übrigens "hello.cc", die Routinen für den OpenCL-Teil sind in "myopencl.cc" und der Kernel heisst "hello_kernel.cc".
Die unterschiedlichen Optimierungsversuche sind in "kernelversionen/" zu finden.
Beim Aufruf kann man mit Option -o die Optimierungsvariante auswählen und mit Option -t noch ein anderes TPB setzen. (Option -? gibt eine kleine Hilfe.)

Fehlersuche (debuggen des OpenCL-Kernels)

Die Suche nach Informationen über OpenCL war etwas schwierig. Einmal sind die Funktionsnamen von OpenCL bei c und c++ unterschiedlich und deshalb schwierig die richtigen zu finden. Dann die Informationsseiten der KhronosGroup sind teilweise extrem unübersichtlich.

Deshalb hier eine kleine Zusammenstellung der wichtigsten Funktionen:

Funktion	Parameter-Beispiele	Beschreibung
int err; err=program.build({device},optionen)	"-cl-nv-verbose" "-cl-nv-maxrregcount=32"	das "verbose" ist nötig wenn man eine Rückmeldung haben will
std::string default_device.getInfo<CL_DEVICE_NAME>()		gibt Name der Grafikkarte als String zurück
std::string program.getBuildInfo<CL_PROGRAM_BUILD_LOG>(device)	default_device	Fehlermeldungen oder Informationen bei erfolgreichem Compilieren
size_t bin_sz; program.getInfo(CL_PROGRAM_BINARY_SIZES,&bin_sz)	&bin_sz	gibt die Grösse des PTX-Codes zurück
char *bin=new char[bin_sz]; program.getInfo(CL_PROGRAM_BINARIES,bin)		gibt den PTX-Code als Text zurück
queue.enqueueNDRangeKernel(kernel,offset,globrange,locrange)	cl::NDRange(xmax), cl::NDRange(TPB)	offset wird nie gebraucht und sollte immer cl::NullRange sein

Einige Funktionen geben einen Fehlercode zurück, 0 wenn erfolgreich, oder eine negative Zahl wenn Fehler. Die Definitionen der Fehlercodes sind in /usr/include/CL/cl.h zu finden. Die Funktion getErrorString(int) in "myopencl.cc" zeigt in etwa diese Definitionen an.

Mit einem "#define __CL_ENABLE_EXCEPTIONS" vor dem "#include <CL/cl.hpp>" hat man die Möglichkeit try-catch-Konstrukte zu verwenden:

  try {
  ...
  } catch(cl::Error &error) {
  ...
  }

  } catch(cl::Error &error) {
  std::cout << error.what() << "(" << error.err() << ") "
            << getErrorString(error.err()) << std::endl;
  std::string name     = default_device.getInfo<CL_DEVICE_NAME>();
  std::string buildlog = program.getBuildInfo<CL_PROGRAM_BUILD_LOG>(default_device);
  std::cerr << "Build log for " << name << ":" << std::endl << buildlog << std::endl;
  return false;
 }

Bei erfolgreichem Compilieren bekommt man mit program.getBuildInfo() in etwa folgenden Ausdruck:

Build log for GeForce GT 625:

ptxas info    : 0 bytes gmem
ptxas info    : Compiling entry function 'hello' for 'sm_21'
ptxas info    : Function properties for hello
ptxas         .     0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas info    : Used 31 registers, 32800 bytes smem, 64 bytes cmem[0]

Nachdem das Compilieren erfolgreich war, kann man beim laufen lassen immer noch Fehler bekommen. Eine typische Fehlermeldung sieht dann so aus:

Fehler10: Kernelaufruf misslungen:
  clFinish(-36) CL_INVALID_COMMAND_QUEUE
device_speicher_freigeben()
terminate called after throwing an instance of 'cl::Error'
  what():  clFinish
Abgebrochen (Speicherabzug geschrieben)

Erzeugen einer PTX-Datei

 size_t bin_sz;
 program.getInfo(CL_PROGRAM_BINARY_SIZES,&bin_sz);
 char *bin=new char[bin_sz];
 program.getInfo(CL_PROGRAM_BINARIES,&bin);
 FILE *fp=fopen("tmp.ptx","w");
 fprintf(fp,"%s\n",bin);
 fclose(fp);
 delete[] bin;

Makefile

Wenn schon defaultmässig eine optimierte Variante verwendet werden soll, kann man diese von kernelversionen/ ins Hauptverzeichnis kopieren und als "hello_kernel.cc" speichern.
Damit das Hauptprogramm weiss welche Optimierung verwendet werden soll, ist jeweils ein entsprechendes #define eingebaut:

#ifdef NUR_OPTIONEN
OPTIMIERUNG 5
#else
__kernel ...
...
}
#endif

Auszug aus dem makefile:

KERNEL = hello_kernel.cc
KERNEL_OPT = kernelversionen/hello_kernel-opt5.cc

all: hello
#all: hello tmp.o

hello: hello.cc hello.h myopencl.cc myopencl.h

check: $(KERNEL_OPT)
        gcc -c -I/usr/local/cuda/include/ -DCHECK $(KERNEL_OPT) -o tmp.o
tmp.o: $(KERNEL)
        gcc -c -I/usr/local/cuda/include/ -DCHECK $(KERNEL) -o tmp.o

Damit der gcc bei den OpenCL-spezifischen Sachen nicht reklamiert, sollte noch die Datei "kernelcheck.h" im Kernel eingefügt werden, aber nur wenn CHECK gesetzt ist:

#ifdef CHECK
#include "kernelcheck.h"
#endif

vollständiges Hello-World mit OpenCL

Ich habe dazu das oben erklärte Beispiel verwendet, von dem ich glaube dass es so einfach wie möglich, aber nicht einfacher ist.

Downloads:

beispiel3_opencl.tar.gz   Hauptprogramm mit Optionen zum alle Optimierungsvarianten auszuprobieren

Fortsetzungen: