[X] OpenMP

[Ben04]

Inhaltsverzeichnis

• 1 OpenMP - Was ist das?

• 2 Compiler Unterstützung

• 3 Syntax

• 3.1 Parallel

• 3.2 Thread ID, single und master

• 3.3 Sections

• 3.4 Barrier

• 3.5 For

• 3.4.1 Verschachtelte Schleifen mit und ohne Collapse

• 3.4.2 Partielle sequentielle Abarbeitung mit ordered

• 3.6 Kurzschreibweise

• 4 Gemeinsamer Speicher

• 4.1 Ciritcal

• 4.2 Flush

• 4.3 Reduction

• 4.4 Locks/Mutex

• 5 Zeitmessung

• 6 Suchen in einem ungeordneten Array

• 7 Akkumulation

• 7.1 STL

• 8 Gibt es noch mehr?

• 9 Quellen

1 OpenMP - Was ist das?

OpenMP ist eine standardisierte Schnittstelle für Fortran, C und C++, welche die Entwickelung von Programmen, welche mehrere Prozessor-Kerne nutzen, vereinfachen soll. Ich werde in diesem Artikel aber nur auf die C++ Schnittstelle eingehen.

2 Compiler Unterstützung

Es gibt mehrere Versionen von OpenMP. Dieser Artikel beschäftigt sich hauptsächlich mit OpenMP 2.5. Ich werde allerdings auch auf ein paar Neuerungen in der Version 3.0 eingehen. OpenMP wird von einer Reihe von Compilern unterstützt, jedoch meistens nicht in den Basisausstattungen dieser. Mit dem GCC gibt es aber mindestens ein Compiler der jedem Leser zugänglich ist.

Folgende Tabelle liefert einen kleinen Überblick.

• GCC >= 4.2 (mit dem -fopenmp Parameter)
• MinGW Technology-Preview-Version (mit dem -fopenmp Parameter)
• VC 2005 Professional und Team System (mit dem /openmp Parameter, Standard und Express Versionen gehen nicht)

Auf der OpenMP-Webseite findet man eine detailliertere Tabelle.

3 Syntax

OpenMP besteht zu einem Großteil aus Compiler-Erweiterungen. Es gibt eine Reihe von Pragma-Anweisungen welche die OpenMP-Funktionalität bereitstellen. Jede von ihnen hat folgenden Syntax:

#pragma omp Anweisung
C++ Code auf den sich die Anweisung bezieht

Neben diesen pragma-Konstrukten gibt es noch eine Reihe von normalen Funktionen. Diese werden im Header omp.h definiert. Alle haben den Präfix omp_.

Jeder Compiler der sich an den C++-Standard hält muss unbekannte Pragma-Anweisungen ignorieren und höchstens mit einer Warnung versehenen. Fehler sich nicht standard-konform. Durch das Benutzen eines Dummy-Headers, ist es deswegen in vielen Fällen möglich eine OpenMP-Anwendung sequentiell auf nicht OpenMP-fähigen Compilern zu übersetzen. Dies ist auch empfehlenswert wenn der Compiler zwar OpenMP unterstützt, allerdings der Zielprozessor keine echte Parallelisierung unterstützt.

In diesem Artikel werde ich versuchen die Anweisungen welche kein sequentielles Fallback-Konstrukt haben zu ignorieren. Sie bieten auch eher selten Funktionalität an, welche man nicht auch anders erhalten kann. Daher ist dies an sich kein großer Verlust.

3.1 Parallel

Die grundlegend Anweisung, ohne die bei OpenMP nichts auskommt, ist die parallel-Anweisung. Sie erzeugt eine Reihe von Threads und führt den nachfolgenden C++ Block mit jedem Thread aus. Dies bedeutet dass der gleiche Code mehrmals gleichzeitig ausgeführt werden kann und zwar von jedem Thread einmal. Wie viele Threads erstellt werden, ist der Implementierung von OpenMP freigestellt. Meistens werden aber so viele Threads erzeugt wie es im Rechner Kerne gibt.

Es ist also in der Regel nicht möglich vorher zu sagen wie oft ein parallel-Block ausgeführt wird.

int main(){
	printf("Hier ist noch alles einfach und sequentiell.\n");
	#pragma omp parallel
	{
		printf("Hallo ");
		printf("Welt!\n");
	}
	printf("Alle Threads haben Hallo gesagt.");
}

Im parallel-Block ist ganz normaler C++-Code erlaubt. Man kann C++-Objekte anlegen und auch Funktionen aufrufen. Sogar Ausnahmen sind erlaubt, solange sie nicht aus dem parallel-Block rausfliegen.

Solange sich die Threads nicht in die Quere kommen, ist das auch schon alles was man machen muss um sämtliche Kerne seines Rechners zu benutzen. Allerdings ist dies leider nur sehr selten der Fall ist. Bereits in diesem einfachen ist es möglich, dass Threads sich streiten. Wenn zwei Threads gleichzeitig etwas ausgeben, ist es nämlich nicht klar welcher zuerst dran kommt. Alles von Funktionieren-wie-erwartet bis hin zu Abstürzen ist eine möglich Konsequenz.

Was das Problem noch komplizierter macht, ist dass das Verhalten nicht vollständig deterministisch ist. Wie viele Threads erzeugt werden kann von Lauf zu Lauf unterschiedlich sein und welcher Thread wann dran kommt kann man noch schlechter vorhersagen. Dies führt dazu, dass manche Bugs nur mit einer gewissen Wahrscheinlichkeit eintreten. Das Debuggen solcher Fehler ist sehr schwer.

An sich müsste ich hier sicherstellen, dass zu jedem Zeitpunkt immer höchstens nur ein Thread etwas auf die Konsole schreibt. Bei printf kann man aber normalerweise davon ausgehen, dass legendlich die Reihenfolge der Ausgaben etwas durcheinander geraten kann. Abstürze gehören eher in den Bereich der Theorie. printf leidet normalerweise weniger unter diesem Problem als std::cout, deswegen werde ich in diesem Artikel ausnahmsweise printf bevorzugen. Im Folgenden werde ich davon ausgehen, dass Ausgaben mit printf einfach funktionieren. Alles andere würde die Beispiele zu kompliziert machen.

Da man nicht wirklich weiß wie oft ein parallel-Block ausgeführt wird, ist eine Parallelisierung recht umständlich ohne weitere OpenMP-Konstrukte. Er eignet sich aber in Regel vorzüglich um threadspezifische Variablen zu initialisieren.

Mit Hilfe der num_threads(n)-Anweisung kann man festlegen wie viele Threads erstellt werden. n muss eine Ganzzahl sein und braucht nicht zum Zeitpunkt der Übersetzung bekannt zu sein.

int main(){
	printf("Wie viele Threads sollen es sein?\n");
	int n;
	scanf("%d", &n);
	#pragma omp parallel num_threads(n)
	{
		printf("Hallo ");
		printf("Welt!\n");
	}
	printf("Alle Threads haben Hallo gesagt.");
}

Wenn man für n eine große Zahl eingibt, kann man mit großer Wahrscheinlichkeit beobachten wie sich die Threads in die Quere kommen. Es folgt eine Beispielausgabe, so wie ich sie beobachten konnte.

Hallo Hallo Welt!
Hallo Welt!
Welt!
Hallo Welt!

Während der eine Thread noch am Schreiben ist, funkt ein anderer ihm einfach dazwischen. Versuchen Sie nach zu vollziehen wie diese Ausgabe zustande kam. Es ist sehr wichtig dieses Problem zu verstehen, denn sonst werden Sie mit OpenMP oder auch Multithread-Programmierung im Allgemeinen nur Ärger haben.

omp_get_max_threads gibt an wie viele Threads im nächsten parallel-Block erzeugt werden.

Bei den meisten OpenMP-Implementierungen werden die Threads nur einmal erstellt. Dies gilt auch wenn es mehrere parallel-Blocks gibt. Zwischen den Blocks gehen die zusätzlichen Threads in einen Wartezustand. Sie werden erst am Ende des Programms zerstört. Die in Regel recht hohen Threaderstellungskosten muss man also nur genau einmal zahlen.

parallel-Blöcke dürfen verschachtelt werden. In der Standardeinstellung erzeugt der verschachtelte Block aber keine neuen Threads.

3.2 Thread ID, single und master

Jeder Thread besitzt eine eindeutige Identifikationsnummer. Es handelt sich dabei um fortlaufende Nummern und nicht etwa um gecastete Zeiger oder ähnliches. Man kann die ID also als Index in einem Array benutzen. Den Thread mit der ID 0 nennt man auch noch Master.

Mit Hilfe von omp_get_thread_num kann man die ID eines Threads bestimmen und durch einen Aufruf von omp_get_num_threads kann man die Gesamtanzahl der Threads, welche vom parallel-Block erzeugt wurden, heraus finden.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());	
		if(omp_get_thread_num() == 0)
			printf("Es gibt %d Threads.\n", omp_get_num_threads());		
	}

}

Wie man an diesem Beispiel sieht, ist es möglich auch in einem parallel-Block, Instruktionen nur einmal aus zu führen. Allerdings gibt es dafür bessere OpenMP-Konstrukte.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());	

		#pragma omp master
		{
			printf("Es gibt %d Threads.\n", omp_get_num_threads());	
		}

		#pragma omp single
		{
			printf("Thread %d war hier.\n", omp_get_thread_num());	
		}

		#pragma omp single
		{
			printf("Danach war Thread %d da.\n", omp_get_thread_num());	
		}
	}
}

Ein master-Block wird nur vom Master-Thread ausgeführt. Der Code ist äquivalent zur if-Abfrage aus dem vorherigen Beispiel. Ein single-Block wird genau einmal ausgeführt, allerdings ist es egal von welchem Thread.

Sämtliche Threads warten am Ende eines single-Blocks bis auch der letzte es bis dahin geschafft hat. Dadurch ist sichergestellt, dass der zweite single Block nicht ausgeführt wird ehe der erste vollständig abgearbeitet wurde. Durch eine nowait Anweisung kann man das Warten unterdrücken.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			printf("Manchmal bin ich als erstes da.\n");	
		}

		#pragma omp single nowait
		{
			printf("Hin und wieder gewinne aber auch ich das Rennen.\n");	
		}
	}
}

Wenn man nowait verwendet dann wird nur noch garantiert, dass die Ausführung des ersten Blocks vor der des zweiten beginnt. Es ist aber gut möglich, dass der zweite vor dem ersten beendet wird. Es ist auch möglich, dass die Ausführung des ersten Threads noch vor der ersten Anweisung im Block unterbrochen wird und es dadurch erscheint, als wenn der zweite Block vor dem ersten gestartet werden würde. In der Praxis kann man also nicht darauf bauen, dass der erste vor dem zweiten gestartet wird, auch wenn dies mit einer hohen Wahrscheinlichkeit so sein wird.

Am Ende eines Master-Block wird nicht gewartet.

Bei verschachtelten parallel-Blöcken werden die IDs durch den innersten Block neu verteilt. In der Standardeinstellung erhalten alle Threads die ID 0 und die Anzahl der Threads wird als 1 gemeldet. Die Idee dahinter ist, dass sich alles auf den innersten Block bezieht. Dieser hat aber keinen Thread erzeugt und folglich gibt es nur den der den Block angestoßen hat. Eine weitere Konsequenz ist, dass single- und master-Blöcke mehrfach ausgeführt werden. Ein verschachtelter parallel-Block wird demnach also in der Standardeinstellung sequentiell abgearbeitet. Mit num_threads(n) kann das erzeugen weiterer Threads erzwingen.

Die Problematik wird durch folgendes Beispiel gut illustriert.

int main(){
	#pragma omp parallel
	{
		#pragma omp parallel
		{
			printf("omp_get_thread_num() = %d,"
				" omp_get_num_threads() = %d\n", 
				omp_get_thread_num(), 
				omp_get_num_threads());	
		}
	}
}

3.3 Sections

Alternativ zu den single-Blöcken kann man auch sections verwenden.

int main(){
	#pragma omp parallel
	{
		#pragma omp sections
		{
			{
				printf("Thread %d war hier.\n", omp_get_thread_num());			
			}
			#pragma omp section
			{
				printf("Aber hier war Thread %d.\n", omp_get_thread_num());		
			}
			#pragma omp section
			{
				printf("und eine weitere section");		
			}
		}
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

Jede der drei Sektionen wird genau einmal ausgeführt. Am Ende des sections-Block warten die Threads ehe sie fortfahren. Auch hier kann man das Warten durch nowait unterdrücken.

Man könnte obiges Beispiel auch mit single ausdrücken.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			printf("Thread %d war hier.\n", omp_get_thread_num());			
		}
		#pragma omp single nowait
		{
			printf("Aber hier war Thread %d.\n", omp_get_thread_num());		
		}
		#pragma omp single nowait
		{
			printf("und eine weitere section");		
		}
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

Bei der single-Variante ist der erste Block mit hoher Wahrscheinlichkeit der erste der gestartet wird. Bei der sections-Variante sind alle Blöcke gleich wahrscheinlich.

3.4 Barrier

Wenn an einer beliebigen Stelle gewartet werden soll so kann man eine barrier-Anweisung benutzen.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());
		#pragma omp barrier
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

In diesem Beispiel ist sicher gestellt, dass alle Threads "Hallo" sagen ehe sich der erste verabschiedet.

3.5 For

Mit Hilfe von OpenMP kann man auch einfache Schleifen parallelisieren.

int main(){
	#pragma omp parallel
	{
		#pragma omp for
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);
	}
}

Anders als bei einer normalen Schleife, werden die Durchläufe hier nicht unbedingt nacheinander abgearbeitet, sondern möglicherweise gleichzeitig und auch in zufälliger Reihenfolge. Dies führt dazu, dass die Ausgabe des obigen Codes nicht immer die gleiche ist. Die Reihenfolge der Zeilen kann beliebig permutiert sein.

Um dies zu bewerkstelligen werden die Durchläufe möglichst gleichmäßig zwischen den Threads verteilt. Zum Beispiel bei zwei Threads, könnte der erste die Durchläufe mit 0<=i<=4 abarbeiten und der andere die restlichen mit 5<=i<=9. Die genaue Verteilung ist dem Compiler frei gestellt. Sie hängt jedoch nur von der Anzahl der Threads und der Anzahl der Schleifendurchläufen ab, das heißt wenn aus irgendeinem Grund ein Thread massive Verspätung hat, dann wird nicht kurzerhand die Arbeit umgeschichtet.

int main(){
	#pragma omp parallel
	{
		// Sorgt für massive Verspätung eines Threads
		#pragma omp single nowait
		{
			for(int i=0; i<INT_MAX; ++i)
				for(int j=0; j<INT_MAX; ++j)
					{}
		}

		// Tja hier warten dann die anderen...
		#pragma omp for
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);
	}
}

Obwohl der eine Thread sehr lange mit dem single-Block beschäftigt ist, kommt keiner der wartenden Threads auf die Idee bei der for-Schleife Überstunden zu leisten. Dies ist nur mit einer dynamischen Arbeitsverteilung möglich.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			for(int i=0; i<INT_MAX; ++i)
				for(int j=0; j<INT_MAX; ++j)
					{}
		}

		#pragma omp for schedule(dynamic)
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);

		// Hier wird erst gewartet wenn die Schleife vollständig abgearbeitet wurde.
	}
}

Wenn ein Thread fertig ist, schaut er ob es noch Arbeit gibt und übernimmt diese dann auch gleich. Ein ähnliches Problem ergibt sich wenn die Durchläufe unterschiedlich lang dauern.

Am Ende eines for-Blocks wird gewartet. Dies kann man wie üblich mit nowait unterdrücken.

Da schedule(dynamic) einen größere Verwaltungskosten als die normale Version hat, kann es sein, dass die normale Version in manchen Fällen schneller ist, obwohl nicht alle Threads immer beschäftigt sind. Die normale Version kommt nämlich ganz ohne Kontakt zwischen den einzelnen Threads zurecht und ist daher sehr schnell. Ich werde weiter unten auf konkrete Implementierungsmöglichkeiten eingehen.

Damit eine Schleife überhaupt parallelisiert werden kann muss sie eine Reihe von sehr strikten Kriterien erfüllen.

• Die Schleifenvariable muss eine Ganzzahl sein. In OpenMP 3.0 sind auch Random-Access-Iteratoren und Zeiger erlaubt.
• Die Abbruchbedingung muss ein Ordnungsvergleich (<,>, <= oder sein und auf der linken Seite muss sich die Schleifenvariable und nur die befinden. Die rechte Seite darf nicht während der Ausführung der Schleife verändert werden.
• Die Inkrementierung muss in immer gleich großen Schritten stattfinden. Dies bedeutet, dass man nur --, ++, += oder -= verwenden darf um die Schleifenvariable zu verändern. Die Größe der Schritte darf nicht während der Ausführung der Schleife verändert werden.
• break und return sind verboten. continue ist erlaubt. goto darf nur verwendet werden wenn man nicht aus der Schleife raus springt.

Folgende Schleifen lassen sich nicht mit OpenMP parallelisieren.

int n = foo();

// Auf der einen Seite der Bedingung muss ein einfaches i stehen.
for(int i=0; i*i<n; ++i)
	printf("%d", i);

// Die Inkrementierungsschritte sind nicht gleichgroß
for(int i=0; i<n; i*=3)
	printf("%d", i);

// Dies ist erlaubt, allerdings ist unspezifiziert wie oft bar und foobar aufgerufen werden.
for(int i=-10; i<n*n*bar(); i+=foobar(n))
	printf("%d", i);

// Man darf nicht aus der Schleife raus springen.
int a[10] = { ... }
for(int j=0; j<10; ++j)
	if(a[j] == n)
		return i;

// n darf nicht in der Scheife verändert werden.
for(int i=0; i<n*10; i+=n)
	++n;

3.4.1 Verschachtelte Schleifen mit und ohne Collapse

OpenMP-for-Schleifen dürfen nicht verschachtelt werden. Man kann allerdings mit einem zusätzlichen parallel-Block Abhilfe schaffen. Dies ist jedoch selten sinnvoll da die Schleife zwar eine OpenMP-Schleife ist, aber dennoch sequentiell abgearbeitet wird.

#pragma omp parallel
{
	#pragma omp for
	for(int i=0; i<10; ++i)
	{
		#pragma omp parallel
		{
			#pragma omp for
			for(int j=0; j<10; ++j)
			{
				printf("%d*%d = %d\n", i, j, i*j);
			}
		}
	}
}

Es wird aber nur die äußere Schleife parallelisiert und die innere sorgt nur für Probleme und Kosten. Besser ist:

#pragma omp parallel
{
	#pragma omp for
	for(int i=0; i<10; ++i)
	{
		// Dieses for ist ein ganz normales C/C++ for
		for(int j=0; j<10; ++j)
		{
			printf("%d*%d = %d\n", i, j, i*j);
		}

	}
}

In OpenMP 3.0 gibt es eine spezielle collapse-Schleife um über ein mehrdimensionales Array zu iterieren. Dies entspricht einem der häufigsten Einatzgebiet von mehrfach verschachtelten Schleifen. Die Schleifen müssen allerdings sehr einfach gehalten sein, damit eine Parallelisierung überhaupt erfolgen kann. Alle Schleifenvariablen müssen völlig unabhängig von einander sein. Es darf nur die innerste Schleife beliebigen Code enthalten. Alle anderen dürfen nur verschachtelte Schleifen enthalten. Das Argument der collapse-Anweisung gibt an wie tief die Schleifen verschachtelt sind.

Dies sind sehr starke Einschränkungen und deswegen deckt das folgende Beispiel an sich auch quasi jede erlaubte Variante einer collapse-Schleife ab.

int main(){
	#pragma omp parallel
	{
		int foo[10][3][7];
		#pragma omp for collapse(3)
		for(int i=0; i<10; ++i)
			for(int j=0; j<3; ++j)
				for(int k=0; k<7; ++k){
					// Nur hier darf beliebiger Code stehen!
					foo[j][k] = 0;
				}

}

3.4.2 Partielle sequentielle Abarbeitung mit ordered

Man kann in eine OpenMP-[i]for* Schleife einen ordered-Block packen. Dieser sorgt dafür, dass ein Teil der Schleife sequentiell abgearbeitet wird.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<=10; ++i){
			double v = i;
			for(int j=1; j<10; ++j)
				v *= i;
			#pragma omp ordered
			{
				printf("%d^10 = %.0f\n", i, v);
			}
		}
	}
}

Die Berechnung der Werte kann in beliebiger Reihenfolge und parallel erfolgen. Die Ausgabe ist jedoch sequentiell, das heißt die Ausgabe ist immer die gleiche.

0^10 = 0
1^10 = 1
2^10 = 1024
3^10 = 59049
4^10 = 1048576
5^10 = 9765625
6^10 = 60466176
7^10 = 282475249
8^10 = 1073741824
9^10 = 3486784401
10^10 = 10000000000

Ein ordered-Block darf nur in einer speziellen ordered-for-Schleife auftreten. Jeder Durchlauf darf nur maximal einem ordered-Block begegnen. Folgendes Beispiel ist demnach falsch.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<10; ++i){
			#pragma omp ordered
			printf("%d^10 =", i);

			double v = i;
			for(int j=1; j<10; ++j)
				v *= i;

			// Hier wird etwas schief gehen da 
			// dies der zweite Block ist.
			#pragma omp ordered
			printf("%.0f\n", v);
		}
	}
}

Dies bedeutet jedoch nicht, dass man nur einen ordered-Block pro Schleife haben darf.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<10; ++i){
			if(i%3 == 0){
				#pragma omp ordered
				printf("%d ist ein Vielfaches von 3\n", i);
			}else{
				#pragma omp ordered
				printf("%d ist kein ein Vielfaches von 3\n", i);
			}
		}
	}
}

In dem Beispiel ist sichergestellt, dass kein Durchlauf zwei verschiedene Blöcke zu sehen bekommt.

3.6 Kurzschreibweise

Wenn ein parallel-Block nur eine for-Schleife enthält, dann gibt es eine Kurzschreibweise.

#pragma omp parallel for
for(int i=0; i<10; ++i)
{
	for(int j=0; j<10; ++j)
	{
		printf("%d*%d = %d\n", i, j, i*j);
	}

}

Analog gibt es auch ein parallel sections-Block.

Folgendes Beispiel berechnet alle Primzahlen kleiner als eine gegebene Zahl.

bool is_prime(unsigned n){
	for(unsigned i=2; i<n; ++i)
		if(n%i == 0)
			return false;
	return true;
}

int main(){
	unsigned n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%u", &n);

	// schedule(dynamic) ist notwendig da is_prime keine
	// konstante Laufzeit besitzt. 
	#pragma omp parallel for schedule(dynamic)
	for(int i=2; i<=n; ++i)
		if(is_prime(i))
			printf("%u ist prim\n", i);
}

4 Gemeinsamer Speicher

Die Zugriffsregeln für Variablen sind die die man erwarten würde. Variablen leben so lange wie der Block in dem sie definiert wurden. Sie können verdeckt werden und man kann auf Variablen einer höheren Ebene zugreifen.

int a;

int main(){
	int b;
	#pragma omp parallel
	{
		int c;
		// a und b werden geteilt, jeder Thread hat sein eigenens c
	}
}

Wenn ein Thread schreibend auf a oder b zugreift dann ergeben sich sämtliche vom Multithreading her bekannte Probleme. OpenMP bietet aber leider nur wenig Revolutionäres um dieses alt bekannte Problem zu lösen.

4.1 Ciritcal

Ein critical-Block ist ein Block in dem sich immer nur genau ein Thread befindet. Wenn bereits ein Thread im Block ist, wartet jeder Thread der hinein will, bis der Block wieder frei wird.

Betrachten Sie folgendes Beispiel:

int main(){
	#pragma omp parallel for
	for(int i=0; i<10; ++i){
		int n;
		#pragma omp critical
		{
			printf("Geben sie eine Zahl ein: ");
			scanf("%d", &n);
		}
		printf("%d * %d = %d\n", i, n, i*n);
	}
}

Hier wird ein critical-Block verwendet um sicher zu stellen, dass der Benutzer nie aufgefordert wird zwei oder mehr Zahlen gleichzeitig einzugeben. Nimmt man das critical weg dann kann es zu folgender Ausgabe kommen.

Geben sie eine Zahl ein: Geben sie eine Zahl ein:

Jeder critical-Block besitzt einen Namen. Wenn zwei Blöcke denselben Namen besitzen, kann nur genau ein Thread in einem der beiden sein. Zwei Blöcke mit verschiedenen Namen beeinflussen sich nicht. Das heißt es ist durchaus möglich, dass 2 Threads in critical-Blöcken sind sofern sie verschiedene Namen haben.

Wenn bei einem critical-Block kein Name angegeben wird, so wie im ersten Beispiel, besteht der Name aus dem leeren Wort. Alle leeren Wörter werden als gleich angesehen. Im Klartext bedeutet dies, dass keine zwei Threads gleichzeitig in irgendeinem namenlosen Block sein können.

Leider sind die Namen global und liegen in keinem Namensraum. Dadurch sind Namenskonflikte möglich. Es ist allerdings sehr unwahrscheinlich, dass dies andere Auswirkungen hat als einfach nur ein langsameres Programm, da legendlich ein Thread warten muss obwohl dies nicht nötig gewesen wäre. Es ist sichergestellt, dass sie normalen Symbolen wie Funktionsnamen nicht in die Quere kommen.

int main(){
	FILE*file = fopen("out.txt", "w");
	#pragma omp parallel for
	for(int i=0; i<10; ++i){
		int n;

		#pragma omp critical(console)
		{
			printf("Geben sie eine Zahl ein : ");
			scanf("%d", &n);
		}

		#pragma omp critical(console)
		{
			printf("%d * %d = %d\n", n, n, n*n);
		}

		#pragma omp critical(file)
		{
			fprintf(file, "%d * %d = %d\n", n, n, n*n);
			fprintf(file, "%d + %d = %d\n", n, n, n+n);
		}

		#pragma omp critical(console)
		{
			printf("%d + %d = %d\n", n, n, n+n);
		}
	}
	fclose(file);
}

In diesem Beispiel ist sicher gestellt, dass nichts auf die Konsole geschrieben wird, falls der Benutzer aufgefordert wurde etwas einzugeben. In der Datei müssen beide zueinander passenden Rechnungen hinter einander ausgegeben werden. Auf der Konsole ist dies nicht garantiert.

4.2 Flush

Viele Rechner besitzen Register und mehrere Level von Speicher-Caches, um den Zugriff auf Variablen zu beschleunigen. Anstatt direkt mit der Variable im Hauptspeicher zu arbeiten wird diese zuerst in einen schnelleren Cache-Speicher kopiert. Anschließend wird mit der Kopie gerechnet und nach einiger Zeit wird die Variable in den Hauptspeicher zurück kopiert. Dies führt dazu, dass es mehrere Versionen von der gleichen Variable gibt welche nicht immer gleich sind.

In einem sequentiellen Programm stellt dies kein Problem dar, da zu jedem Zeitpunkt bekannt ist wo sich die rezenteste Version der Variable befindet. Jedoch bei mehreren Threads ist dies nicht mehr so einfach. Der eine weiß nicht wann der andere fertig ist und wann er die Werte synchronisieren muss.

Um dieses Problem zu lösen bietet OpenMP die flush-Anweisung an. Sie synchronisiert die lokale Version des Threads mit der des Hauptspeichers. Dies bedeutet, dass jeder Thread nachdem er fertig ist, eine geteilte Variable zu verändern, diese auch flushen muss. Ein Thread der lesend auf eine Variable zugreift will, muss die Variable auch flushen damit er sicher ist die neuste Version zu Gesicht zu bekommen. Wenn ein Thread auch mit einer veralteten Version arbeiten kann, braucht er die Variable nicht zu flushen. Wenn sicher ist, dass eine Variable nicht in einem parallel-Block verändert wird, dann kann man sich das flush auch sparen da es alle Versionen gleich sind, ganz egal wo sie sich befinden.

Folgendes Programm berechnet das Faktoriell einer Zahl.

int main(){
	unsigned n, fact = 1;
	printf("Geben sie eine Zahl ein: ");
	scanf("%u", &n);

	#pragma omp parallel
	{
		unsigned local_fact = 1;

		// Auf n kann man ohne flush zugreifen,
		// da es nicht verändert wird.
		#pragma omp for
		for(int i=1; i<=n; ++i){
			// local_fact braucht nicht geflusht 
			// zu werden, da es nicht geteilt wird.
			local_fact *= i;
		}

		#pragma omp critical
		{
			// Zuerst holt der Thread sich die neuste Version
			// aus dem Hauptspeicher.
			#pragma omp flush(fact)

			// Danach wird die Variable verändert.
			fact *= local_fact;

			// Anschließend wird die Variable wieder in den
			// Hauptspeicher geschrieben, damit andere Threads
			// sich die neuste Version von dort laden können.
			#pragma omp flush(fact)
		}	
	}

	printf("%u! = %u\n", n, fact);
}

Ein fehlendes flush bedeutet nicht, dass die Variablen nicht synchronisiert werden. Es steht der Implementierung frei dies zu tun wann immer sie will. Dies führt leider oft zu nur sehr schwer reproduzierbaren Bugs. Zum Beispiel ist es sehr wahrscheinlich, dass das Betriebssystem die Rechnerzeit anders auf die Threads verteilt, wenn die Hardware unter voller Last läuft, als wenn es nur mit kleinen Testeingaben konfrontiert wird. Es ist daher sinnvoll im Zweifelsfall lieber ein flush zuviel zu schreiben.

4.3 Reduction

Mit reduction lassen sich Berechnung vom Typ s += f(1)+f(2)+...+f(n) parallelisieren. Jeder Thread rechnet ein paar Terme aus und summiert die auf welche er berechnet hat. Nachdem alle Threads fertig sind zählt jeder Thread sein Zwischenergebnis zur Variable s hinzu.

Diese Zwischenergebnis wird in einer Variable mit dem gleichen Namen wie s gespeichert. Diese Hilfsvariable wird im reduction-Block angelegt und mit dem neutralen Element (0 im Fall der Addition) initialisiert. Man sollte sich aber nicht auf diese Verhalten verlassen. Wenn ein Programm sequentiell übersetzt wird, das heißt ohne OpenMP Unterstützung, dann gibt es keine Zwischenergebnisse und deswegen hat s beim ersten Durchlauf auch nicht unbedingt den Wert des neutralen Elements.

Neben + lassen sich auch -, *, &, |, ^, && und || auf diese Weise parallelisieren. Es muss bei den Operatoren jedoch um ihre Buildin-Version handeln.

• +, -, | und ^ besitzen das neutrale Element 0.
• * besitzt das neutrale Element 1.
• & besitzt das neutrale Element ~0. (Bei ~0 sind alle Bits gesetzt sind.)
• && besitzt das neutrale Element true.
• || besitzt das neutrale Element false.

Bei Gleitkommazahlen hängt die Genauigkeit von der Auswertungsreihenfolge ab. Die Reihenfolge ist nicht spezifisiert, also auch nicht die Genauigkeit des Ergebnisses. Es ist nicht einmal garantiert, dass es immer das gleiche Ergebnis sein wird. Gleitkommazahlen in Kombination mit einer Reduktion haben einen gewissen nicht deterministischen Teil im Ergebnis.

In der Fortan Version von OpenMP gibt es zusätzlich noch Min- und Max-Reduktionen. Die wurden aus mir unerklärlichen Gründen in der C/C++ Version vergessen.

Folgendes Beispiel bestimmt ob eine Zahl prim ist.

int main(){
	int n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%d", &n);

	// Der Initialwert von is_prime ist wichtig da dieser auch
	// in die Berechnung einbezogen wird.
	bool is_prime = true;

	#pragma omp parallel reduction(&&: is_prime)
	{
		// jeder bekommt sein eigenes is_prime,
		// welches mit true, dem neutralen Element von &&,
		// initialisiert wurde.

		#pragma omp for
		for(int i=2; i<n; ++i)
			if(n%i == 0)
				is_prime = false;

		// Alle is_prime (auch das aus der main Funktion) werden 
		// mittels && verknüpft. Im Klartext bedeutet dies, dass 
		// falls eines false ist, ist das Resultat false.
	}

	if(is_prime)
		printf("%d ist prim\n", n);
	else
		printf("%d ist nicht prim\n", n);

}

Am Ende des parallel Blocks wird folgende Rechnung ausgeführt:

is_prime = is_prime && thread[0].is_prime && thread[1].is_prime && ...  && thread[thread_num-1].is_prime

Bei der thread[j]-Notation handelt es sich legendlich um Pseudocode.

Auch hier gibt es eine Kurzschreibweise. Folgendes Beispiel berechnet das Faktoriell einer Zahl.

int main(){
	unsigned n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%u", &n);

	unsigned fact = 1;

	#pragma omp parallel for reduction(*: fact)
	for(int i=2; i<n; ++i)
		fact *= i;

	printf("%u! = %u\n", n, fact);
}

reduction erlaubt es dem Compiler eine Reihe von Optimierungen durchzuführen, welche bei der äquivalenten Formulierungen mit critical und flush nur sehr schwer möglich sind.

4.4 Locks/Mutex

OpenMP bietet auch Mutexe an. Hier weicht OpenMP an sich nicht von allgemein bekanten Konzepten ab. Es gibt einen Lock-Type namens omp_lock_t. Jeder Mutex ist entweder offen oder geschlossen. Folgende Funktionen werden bereitgestellt um damit zu arbeiten:

• void omp_init_lock(omp_lock_t*)
  Initialisiert den Mutex. Der Mutex ist danach offen.
• void omp_destroy_lock(omp_lock_t*)
  Zerstört den Mutex.
• void omp_set_lock(omp_lock_t*)
  Versucht den Mutex zu schließen. Wurde der Mutex bereits von einem anderen Thread geschlossen so wird gewartet bis, dass dieser ihn wieder geöffnet hat.
• void omp_unset_lock(omp_lock_t*)
  Öffnet den Mutex wieder
• int omp_test_lock(omp_lock_t*)
  Versucht den Mutex zu schließen. Wurde der Mutex bereits von einem anderen Thread geschlossen so wird 0 zurück gegeben, ansonsten 1.

omp_lock_t ist ein nicht rekursives Mutex. Das heißt ein Thread darf einen Mutex nicht mehrmals schließen. Er würde in dem Fall nämlich auf sich selbst warten. Es gibt aber auch eine rekursive Variante und die heißt omp_nest_lock_t. Die Funktionen um damit umzugehen heißen ähnlich und funktionieren analog.

• void omp_init_nest_lock(omp_next_lock_t*)
• void omp_destroy_nest_lock(omp_next_lock_t*)
• void omp_set_nest_lock(omp_next_lock_t*)
• void omp_unset_nest_lock(omp_next_lock_t*)
• int omp_test_nest_lock(omp_next_lock_t*)

Da diese Schnittstelle alles anderes als modernes C++ ist, bietet es sich an diese Funktionen zu wrappen.

class Mutex
{
public:
	Mutex()				{omp_init_lock(&lock);}
	~Mutex()			{omp_destroy_lock(&lock);}
	void lock()			{omp_set_lock(&lock);}
	void unlock()			{omp_unset_lock(&lock);}
	bool try_to_lock()		{return omp_test_lock(&lock);}
private:
	Mutex(const Mutex&);
	Mutex&operator=(const Mutex&);
	mp_lock_t lock;
};

class ScopedLock
{
public:
	ScopedLock(Mutex&lock)		:lock(lock){lock.lock();}
	~ScopedLock()			{lock.unlock();}
private:
	ScopedLock(const ScopedLock&);
	ScopedLock&operator=(const ScopedLock&);
	Mutex&lock;
};

5 Zeitmessung

Neben den Funktionen welche direkt mit der Multithread-Programmierung bietet OpenMP noch zwei nützliche Funktionen zur Zeitmessung an.

• double omp_get_wtime();
  Gibt die momentane Zeit in Sekunden zurück. Es ist nicht definiert wann der Zeitpunkt 0 ist. Es ist legendlich bekannt, dass es nicht während dem Lauf eines Programms verändert wird und in der Vergangenheit liegt. Dies ist aber nicht weiter schlimm, wenn man die Funktion nur für Zeitmessungen einsetzen will. Die Genauigkeit dieser Funktion ist nicht spezifisiert und darf leicht von Thread zu Thread abweichen. Man kann die Genauigkeit aber mit omp_get_wtick in Erfahrung bringen.
• double omp_get_wtick();
  Gibt den kleinstmöglichen Abstand zwischen 2 verschiedenen von omp_get_wtime zurückgegebenen Zeitpunkten an. Der Rückgabewert von omp_get_wtick darf von Thread zu Thread variieren.

6 Suchen in einem ungeordneten Array

Das Suchen in einem unsortierten Array lässt sich sehr leicht mit OpenMP beschleunigen.

// Ohne OpenMP
int find_first(int arr[], int size, int value){
	for(int j=0; j<size; ++j)
		if(arr[j] == value)
			return j;
	return -1;
}

// Mit OpenMP
int find_any(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel for
	for(int j=0; j<size; ++j)
		if(arr[j] == value){
			i = j;
			#pragma omp flush(i)		
		}
	return i;
}

Das aktuelle i braucht nicht zuerst per flush-Anweisung geladen zu werden, da der Wert eh nur überschrieben werden würde.

Wie man bereits am Namen erkennt, ist bei der parallelen Version nicht definiert welchen Wert gefunden wird, wenn es mehrere gleiche Werte im Array gibt. Wenn dies jedoch nötig ist, kann man dieses Problem leicht mit ordered beheben.

int find_first(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel for ordered
	for(int j=size-1; j>=0; --j)
		if(arr[j] == value){
			#pragma omp ordered
			i = j;
			#pragma omp flush(i)		
		}
	return i;
}

Eines der Nachteile der parallelen Version ist, dass immer alle Elemente angeschaut werden. Dies ist aber nicht immer notwendig denn man kann die Schleife abbrechen sobald ein Element gefunden wurde. OpenMP bietet allerdings keine Möglichkeit dies zu tun. Man kann versuchen die Schleife nach zu bauen.

int find_any(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel
	{
		const int begin = size*omp_get_thread_num() / omp_get_num_threads();
		const int end = size*(omp_get_thread_num()+1) / omp_get_num_threads();

		for(int j=begin; j<end; ++j){
			if(arr[j] == value)
				i = j;
			#pragma omp flush(i)
			if(i != -1)
				break;
		}
	}
	return i;
}

Dieser Code dürfte recht nahe an dem sein, was der Compiler aus einer OpenMP-Schleife mit statischer Verteilung macht. In der Regel sollte man aber die OpenMP-for-Schleife bevorzugen. Diese ist weit weniger fehleranfällig. So kann es beim Beispiel oben zum Beispiel zu Überläufen bei der Berechnung von begin und end kommen.

Ein weiteres Probleme ist, dass jeder Durchlauf den Wert von i laden muss. Dies kann je nach Architektur sehr langsam sein und somit den potentiellen Geschwindigkeitsgewinn wieder weg machen. Um dies zweifelsfrei fest zu stellen muss man die Zeit messen.

7 Akkumulation

Die Akkumulation ist die Standardanwendung der reduction-Anweisung. Es geht darum alle Elemente einer Liste auf irgendeine Weise zu kombinieren. Das Addieren aller Zahlen einer Liste ist ein Beispiel welches mit einer solchen Anweisung gelöst werden kann.

int accumulate(int arr[], int size){
	int sum = 0;
	#pragma omp parallel for reduction(+:sum)
	for(int i=0; i<size; ++i)
		sum += arr;
	return sum;
}

Leider handelt es sich jedoch nicht immer um einen Buildinoperator. Sei [i]Matrix* eine Klasse welche eine 2-dimensionale Matrix beschreibt und für welche der *= Operator entsprechend überladen wurden. Der Defaultkonstruktor initialisiert die Matrix auf die Einheitsmatrix.

class Matrix{
private:
	int 
		a,b,
		c,d
	;
public:
	Matrix():a(1), b(0), c(0), d(1){}
	Matrix&operator*=(const Matrix&right){
		Matrix left = *this;
		a = left.a*right.a + left.b*right.c;
		b = left.a*right.b + left.b*right.d;
		c = left.c*right.a + left.d*right.c;
		d = left.c*right.b + left.d*right.d;
		return *this;
	}
};

Matrix multiply(Matrix arr[], int size){
	const unsigned thread_count = omp_get_max_threads();
	Matrix*sub_prod = new Matrix[thread_count]; 
	#pragma omp parallel num_threads(thread_count)
	{
		const unsigned id = omp_get_thread_num();
		const int begin = size*id/thread_count;
		const int end = size*(id+1)/thread_count;

		for(int j=begin; j<end; ++j)
			sub_prod[id] *= arr[j];
	}
	Matrix prod;
	for(int j=0; j<thread_count; ++j)
		prod *= sub_prod[j];
	delete[]sub_prod;
	return prod;
}

Die Idee ist die Assoziativität der Multiplikation auszunutzen um das Produkt in mehrere Teilprodukte zu zerlegen und diese zuerst von jedem Thread einzeln auswerten zu lassen.

Zuerst wird bestimmt wie viele Threads überhaupt zur Verfügung stehen. Anschießend wird ein Array angelegt welches einen Eintrag für jeden Thread besitzt. In diesem Array werden die Zwischenresultate abgelegt. Die Matrizenmultiplikation ist nicht kommutativ, daher verwende ich zur Berechnung der Teilprodukte keine OpenMP-for-Schleife. Diese garantiert nämliche keine bestimmte Reihenfolge. Eine ordered OpenMP-for-Schleife ist zu steif da nicht alle Multiplikation sequentiell ablaufen müssen. Nachdem alle Teilprodukte berechnet wurden werden sie multipliziert und das Ergebnis wird zurückgegeben.

Ich habe keinen STL-vector verwendet da dieser in der Regel nicht dafür ausgelegt ist, mit mehreren Threads klar zu kommen. In diesem speziellen Fall ist es jedoch höchst unwahrscheinlich, dass dies ein Problem darstellen würde. Man könnte natürlich einen Smartpointer verwenden um die Speicherfreigabe zu garantieren. Allerdings ist das im Fall einer Ausnahme das kleinere Problem.

Die einzige Operationen die eventuell eine Ausnahme werfen könnte, sind die Speicherallokierung selbst und der *= Operator. Wenn die Speicherallokierung fehlschlägt, geht kein Speicher verloren. Der *= Operator unserer Matrixklasse wirft natürlich keine Ausnahmen allerdings ist dies bei einer allgemeinen Matrixklasse nicht mehr selbstverständlich. In diesem Fall muss man allerdings sicher stellen, dass die Ausnahme sicher aus dem parallel-Block hinaus propagiert wird. Dies ist ein wesentlich größeres Problem als das kleine Speicherleck.

7.1 STL

Im Forum kommt es immer wieder zu Diskussionen zum Thema Parallelisierung und wie man die einsetzt ohne Programme von Grund auf neu zu bauen. Oft kommt dann die Idee auf die Algorithmen der STL zu parallelisieren und den Benutzer-Code unverändert zu lassen. Dies ist allerdings schwer um zu setzen. Im folgenden werde ich zeigen wie man std::accumulate parallelisieren kann und welche Probleme es dabei gibt.

Es sei folgende recht einfache Funktion gegeben:

template<class Iter, class T>
T accumulate(Iter begin, Iter end, T sum){
	for(; begin!=end; ++begin)
		sum += *begin;
	return sum;
}

Unser Ziel ist es diese Funktion durch Benutzung von OpenMP zu beschleunigen. Auf den ersten Blick erscheint dies wie eine sehr einfache Aufgabe, dies ist sie aber nicht!

Zuerst müssen wir sicherstellen, dass die zu Grunde liegenden +-Operation parallel ausgeführt werden kann. Speichert der Operator zum Beispiel Zwischenwerte in einem statischen Array, dann können wir eine Parallelisierung ohne die Operation selbst zu verändern bereits vergessen. Es ist auch möglich, dass Werte vorberechnet werden müssen und der Operator auf diese zugreift, das heißt es muss etwas initialisiert werden. Zum Beispiel eine Primzahlliste oder das Pascalsche Dreieck zur Bestimmung großer Binomialkoeffizienten sind gute Kandidaten. Dies muss vor dem ersten Aufruf passiert sein. Es darf nicht erst bis zur ersten Verwendung gewartet werden ehe man Initialisierungsarbeit erledigt.

Die +-Operation keine Ausnahme werfen, da C++ leider keine Möglichkeit bietet diese über Threadgrenzen hinaus zu propagieren.

Des Weiteren benötigt man die Assoziativität der +-Operation. Die Assoziativität ist notwendig damit die Summe in Teilsummen aufgespalten werden kann. Die Kommutativität ist notwendig, da nicht klar ist in welcher Reihe die Threads fertig werden. Dies erscheint offensichtlich, allerdings ist das eine zusätzliche Einschränkung gegenüber der sequentiellen Version. Oft ist sie auch nicht gegeben. Zum Beispiel ist die +-Operation von Gleitkommazahlen nicht assoziativ da die Genauigkeit von der Auswertungsreihenfolge abhängt. Normalerweise sagt, man bei Gleitkommazahlen jedoch, dass die Genauigkeit nicht wichtig ist und nimmt einfach an, dass die Assoziativität gegeben ist.

Ohne die Kommutativität kann man auskommen, wenn es sein muss wie das vorherige Beispiel zeigt. Wenn sie allerdings zur Verfügung steht wird die Implementierung einfacher. Ich gehe im Folgenden davon aus, dass sie gegeben ist.

C++-Templates bieten leider praktisch keine Möglichkeiten diese Bedingungen sicher zu stellen. Die Dokumentation der Funktion ist also die einzige Stelle wo man sie angeben kann. Es bleibt nur zu hoffen, dass der Benutzer der Funktion die Dokumentation auch gründlich ließt und versteht...

Ist sicher gestellt, dass die +-Operation parallel ausgeführt werden kann müssen wir zusehen, dass wir wahlfreien Zugriff auf die Liste von Elementen haben. Dies ist nötig da es einem OpenMP-for frei steht in welcher Reihenfolge er die Schleifenduchläufe ausführen will. Ist dies nicht gegeben, kann man entweder die sequentielle Version benutzen oder die Elemente zu erst kopieren. Diese Kopie muss sequentiell erfolgen, da es keinen wahlfreien Zugriff gibt. Was besser ist hängt wieder von der +-Operation ab.

Bei Ganzzahlen ist die sequentielle Version deutlich überlegen, da die Addition wahrscheinlich so schnell ist wie das Kopieren einer Zahl. Das einfache Ausrechnen ist also in etwa genauso schnell wie das Vorbereiten der parallelen Berechnung. Bei der Matrizenmultiplikation ist es allerdings besser zu kopieren. Werden Matrizen unterschiedlicher Dimensionen multipliziert so lohnt es sich unter Umständen viel mehr, zuerst dynamisch die optimale Auswertungsreihenfolge zu bestimmen als einfach nur zu parallelisieren. Ich werde im Folgenden davon ausgehen, dass es sich lohnt zu kopieren.

Der Code sieht bis hierhin in etwa folgendermaßen aus:

template<class Iter, class T>
T accumulate(Iter begin, Iter end, T init){
	return accumulate_check_category(begin, end, init, std::iterator_traits<Iter>::iterator_category());
}

template<class Iter, class T>
T accumulate_check_category(Iter begin, Iter end, T init, std::random_access_iterator_tag){
	return accumulate_parallel(begin, end-begin);
}

template<class Iter, class T>
T accumulate_check_category(Iter begin, Iter end, T init, std::random_access_iterator_tag){
	typename std::iterator_traits<Iter>::value_type T;
	std::vector<T>buffer(begin, end);
	return accumulate_parallel(buffer.begin(), buffer.size());
}

template<class Iter, class T>
T accumulate_parallel(Iter begin, int len, T init){
	// ... 
}

Die eigentliche Parallelisierung ist ein typischer Reduktionsfall. Die entsprechende OpenMP-Anweisung kann aber nur im Fall eines buildin Typen verwendet werden. Da durch die Parallelisierung Kosten entstehen lohnt es sich nicht unter einer Gewissen Grenze zu parallelisieren. Wie viele Element es geben muss, damit sich das Ganze lohnt hängt von der +-Operation, vom Prozessor und dem verwendeten Compiler ab.

Die OpenMP-Answeisung dürft bei einigen Compilern schneller sein als ihr Nachbau. Beim Nachbau muss man zwischen Nachbau mit Kommutativität und ohne Kommutativität unterscheiden. Bei der OpenMP-Anweisung braucht man dies nicht, da der buildin +-Operator kommutativ ist. Dies macht also mindestens 3 Konstanten, welche man praktisch nur empirisch bestimmen kann, die man in den Code einbauen muss. Da ich die Kommutativität angenommen habe sind es bei mir nur 2. Ich nehme einfach mal 50000 für beide an.

Ein weiteres Problem das sich stellt ist die Frage ob die Laufzeit des Operator+ von seinen Operanten abhängt oder nicht. Ob die for-Anweisung ein schedule(dynamic) benötigt oder nicht hängt davon ab. Will man alle Fälle abdecken so muss man jetzt schon 6 Konstanten bestimmen. Ich gehe mal davon aus, dass es in etwa immer gleich schnell ist. Ein Gegenbeispiel wäre die Multiplikation von Matrizen unterschiedlicher Dimensionen.

template<class Iter, class T>
T accumulate_parallel(Iter begin, int len, T init){
	accumulate_parallel_reduction(begin, len, init, boost::is_arithmetic<typename std::iterator_traits<T>::value_type>::type());
}

template<class Iter, class T>
T accumulate_parallel_reduction(Iter begin, int len, T sum, boost::true_type){
	if(len < 50000){
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}else{
		#pragma omp parallel reduction(+:ret) for
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}
	return sum;
}

template<class Iter, class T>
T accumulate_parallel_reduction(Iter begin, int len, T sum, boost::false_type){
	if(len < 50000){
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}else{
		#pragma omp parallel 
		{
			T sub_sum = T();
			for(int i=0; i<len; ++i)
				sub_sum += *(begin+i);
			#pragma omp critical
			{
				#pragma omp flush(sum)
				sum += sub_sum;
				#pragma omp flush(sum)
			}
		}
	}
	return sum;
}

Ich habe auf meinem 4-Kern-Prozessor und dem GCC-4.2 ein paar Tests gemacht um zu bestimmen wie groß bei mir diese Konstanten liegen. Unter etwa 50000 int-Elementen ist die sequentielle Version bei mir im Vorteil. Bei mehr Elementen ist die parallelisierte Version besser. Erstaunlich ist, dass die Reduktions-Answeisung in etwa genau so schnell ist wie ihr Nachbau. Ich gehe mal davon aus, dass die GCC-Entwickler den OpenMP Code noch nicht wirklich optimiert haben und in dieser Version vor allem auf korrektes Verhalten Wert gelegt haben.

8 Gibt es noch mehr?

In diesem Artikel habe ich nicht alle OpenMP-Anweisungen vorgestellt. Ich habe mich legendlich auf die mir am wichtigsten erscheinenden beschränkt. Wer OpenMP wirklich produktiv einsetzen will, dem empfehle ich die Spezifikation durch zu lesen. Für eine Spezifikation ist diese erstaunlich leicht verständlich und es sind auch viele Beispielanwendungen darin enthalten.

9 Quellen

• http://bisqwit.iki.fi/story/howto/openmp/
• http://www.openmp.org/
• http://algo2.iti.uni-karlsruhe.de/singler/mcstl/

[Ben04]

Sorry für den neuen Thread aber beim alten kriege ich keinen Beitrag mehr gepostet welcher mehr als 2 cpp-Blocks enthält.

==============================

Inhaltlich dürfte sich an der jetzigen Version nicht mehr viel ändern. Ich werde noch vielleicht ein paar Sachen umformulieren.

Ich weiß nicht wie die Unterstützung bei nicht freien Compilern ist. Wenn da jemand ein paar Details kennt, nur her damit.

Andere Kommentare sind natürlich auch erwünscht.

[sothis_]

hab mich noch nie mit openmp beschäftigt, kann daher nichts zum inhalt sagen, ausser das es sehr interessant war und lust darauf gemacht hat mich damit mal auseinander zu setzen (auseinanderzusetzen? :S). von daher

[Artchi]

Der MSVC2005 unterstützt OpenMP 2.0 mit dem /openmp-Parameter. (weiß aber nicht, ob auch in der Express Edition)

[Jester]

Mir fehlt ein bißchen Information, die über eine Dokumentation hinaus geht. Du sagst beispielsweise welche Arten von Schleifen so nicht parallelisiert werden können. Da wäre es schön, wenn Du auch noch ein Sätzchen dazu schreiben könntest wie man sie eben doch parallelisieren kann. (beispielsweise die suche lässt sich ja parallelisieren, indem man die fundstelle einfach speichert und die suche weiterlaufen lässt, die i*=3 schleife lässt sich durch Neuberechnung des Index ja auch parallelisieren usw). Auch eine kleine Anwendung, die über parallele print-statements hinausgeht wäre schön.

merkt OpenMP eigentlich, dass sich so eine Schleife nicht parallelisieren lässt? Oder ist da der Programmierer alleine zuständig?

[Ben04]

Artchi schrieb:

Der MSVC2005 unterstützt OpenMP 2.0 mit dem /openmp-Parameter. (weiß aber nicht, ob auch in der Express Edition)

Danke

Jester schrieb:

beispielsweise die suche lässt sich ja parallelisieren

Meinst du wie folgendermaßen:

int p = -1;
#pragma omp parallel for
for(int i=0; i<n; ++i)
  if(a[i] == x){
    p = i;
    #pragma omp flush(p)
  }
// in p steht jetzt ein Index sodass a[p] == x oder -1 wenn es kein Element gibt.

Wenn ja dann halte ich das für nicht all zu kompliziert wenn man die Idee hinter der Reduktion verstanden hat. Da kann man selber drauf kommen. Desweiteren braucht man ein flush (oder atomic) und um geteilten Speicher geht es erst später. Selbst bei einem Funktor hat man geteilten Speicher. Das Hauptargument war ja, dass man nicht so ohne weiteres rausspringen kann. Das wäre also schon recht weit off Topic.

Jester schrieb:

die i*=3 schleife lässt sich durch Neuberechnung des Index ja auch parallelisieren

Willst du

for(int i=1; i<100; i*=3){...}

als

for(int j=0; j<log(3, 100); ++j){
  int i = exp(3, j);
  ...
}

umschreiben? Das erscheint mir in der Regel als nur wenig sinnvoll da ... schon sehr lange dauern muss damit die parallele Version schneller ist als die sequentielle. Bedenke, dass du wegen des exponentiellen Charakters (bei ints) eh nie viele Durchläufe kriegst.

Das exp ist auch nicht wirklich billig und mir fällt jetzt kein einfacher Ansatz ein das exp weg zu kriegen. Höchstens:

int end = log(3, 100);
#pragma omp parallel
{
  int n = omp_get_num_threads(), t = omp_get_thread_num();
  int my_begin = t*end/n, my_end = (t+1)*end/n;

  int i = exp(3, my_begin);
  // Hier fehlt kein pragma
  for(int j=my_begin; j<my_end; ++j, i *= 3)
    ...
}

Da muss man auf darauf achten, dass es bei (t+1)*end keinen Overflow gibt und dann muss ... noch eine etwa konstante Laufzeit haben.

Außer ich stehe gerade auf dem Schlauch, gibt es keine wirklich einfache und sinnvolle Möglichkeit diese Schleife zu parallelisieren und alles andere hat meiner Meinung nach in einer Beschreibung des Grundkonstrukts OpenMP-for nichts zu suchen. Eventuell könnte man noch ein Kapitel anhängen wie man komplexere Schleifen parallelisiert allerdings habe ich da nicht wirklich ausreichend Erfahrung um zu Entscheiden was da wichtig ist, das heißt rein muss und was nicht.

Den Trick den ich da verwendet habe um die Schleife umzuschreiben könnte ich vielleicht doch einbauen. Man kann ihn auch benutzen um bei der sequentiellen Suche frühzeitig abzubrechen.

Jester schrieb:

Auch eine kleine Anwendung, die über parallele print-statements hinausgeht wäre schön.

Wenn du eine gute Idee für eine Anwendung hast kann ich das machen. Meiner Meinung nach sollte die aber:

• Eine einfache Problemstellung haben da das Hauptgewicht ja auf OpenMP liegt und nicht auf der Anwenung selbst.
• Sinnvoll. Wenn man eine x-beliebige Schleife aus Spaß an der Freude parallelisiert dann kann das gerne nach hinten los gehen, das heißt langsamer als die sequentielle Version sein.

Ich mach lieber ein offensichtlich hinkendes print Beispiel als eine schlechte an den Haaren herbeigezogene Anwendung.

Jester schrieb:

merkt OpenMP eigentlich, dass sich so eine Schleife nicht parallelisieren lässt? Oder ist da der Programmierer alleine zuständig?

Entweder die Schleife wird parallelisiert oder es gibt eine Compilerfehlermeldung.

[Jester]

Ben04 schrieb:

Wenn ja dann halte ich das für nicht all zu kompliziert wenn man die Idee hinter der Reduktion verstanden hat. Da kann man selber drauf kommen. Desweiteren braucht man ein flush (oder atomic) und um geteilten Speicher geht es erst später. Selbst bei einem Funktor hat man geteilten Speicher. Das Hauptargument war ja, dass man nicht so ohne weiteres rausspringen kann. Das wäre also schon recht weit off Topic.

Natürlich ist es nicht unheimlich schwer. Aber es zeigt doch recht gut die Denkweise, die hinter parallelen Algorithmen steht: Lieber mal was zuviel machen (nach auffinden des Elements noch weitersuchen) als dadurch das Skalieren mit wachsender Prozessoranzahl zu verlieren. Das mit dem Speicher ist natürlich ein gewisses Argument... aber ohne geht's wohl schlecht.

Jester schrieb:

die i*=3 schleife lässt sich durch Neuberechnung des Index ja auch parallelisieren

Willst du

for(int i=1; i<100; i*=3){...}

als

for(int j=0; j<log(3, 100); ++j){
  int i = exp(3, j);
  ...
}

umschreiben? Das erscheint mir in der Regel als nur wenig sinnvoll da ... schon sehr lange dauern muss damit die parallele Version schneller ist als die sequentielle. Bedenke, dass du wegen des exponentiellen Charakters (bei ints) eh nie viele Durchläufe kriegst.

Da geht's aber trotzdem um's Prinzip bei der Parallelisierung. Das würde ich nicht hinter "int ist ja klein" verstecken. Nimm's halt mit *=2, dann kannste es mit nem shift noch bequem machen. Es geht (mir zumindest) um die Technik mit der sich möglicherweise sowas eben doch auflösen lässt.

Evtl. könnte man die Schleifenvariable auch über ein tablelookup ermitteln und die lookuptable halbwegs gescheit vorberechnen (quadrieren und die (größerer werdenden) Lücken parallel auffüllen). Aber das ist vielleicht wirklich ein bißchen kompliziert.

Das manuelle Aufsplitten des workload in Teilpakete halte ich für kritisch... das kann man zwar machen (wie du ja sagst, wenn "..." konstante Laufzeit hat), aber eigentlich möchte man das ja nicht tun, sondern das Loadbalancing OpenMP überlassen.

• Eine einfache Problemstellung haben da das Hauptgewicht ja auf OpenMP liegt und nicht auf der Anwenung selbst.
• Sinnvoll. Wenn man eine x-beliebige Schleife aus Spaß an der Freude parallelisiert dann kann das gerne nach hinten los gehen, das heißt langsamer als die sequentielle Version sein.

Akkumulation? Also Summe aller Werte in einem Array? Oder deren Produkt? Allerdings brauchste dafür wieder Speicher... aber ohne kann man halt kaum was interessantes machen.

Zumindest erfüllt es Deine beiden Anforderungen. Außerdem ist es tatsächlich ein absolutes Standardproblem bei parallelen Algorithmen, das man auch wirklich häufig braucht.

[Ben04]

Bin Ich der einzige der Probleme hat hier etwas posten zu können?

sorry für den unkonventionellen post aber anders krieg ich das jetzt nicht hin : http://rafb.net/p/nCtJPJ58.html

Das Forum scheint sich in ne endlos schleife zu verabschieden wenn ich versuche das da zu posten oder in der vorschau an zu sehen. (nicht nur bei dem Post)

[Jester]

sag bitte einfach bescheid wenn du kein feedback haben möchtest. ich lasse es dann. ansonsten verkneif dir bitte deinen imho ziemlich unverschämten ton.

[GPC]

Ben04 schrieb:

Bin Ich der einzige der Probleme hat hier etwas posten zu können?

sorry für den unkonventionellen post aber anders krieg ich das jetzt nicht hin : http://rafb.net/p/nCtJPJ58.html

Das Forum scheint sich in ne endlos schleife zu verabschieden wenn ich versuche das da zu posten oder in der vorschau an zu sehen. (nicht nur bei dem Post)

Ich hab's mal weitergeleitet.

[darthdespotism]

Ben04 schrieb:

Bin Ich der einzige der Probleme hat hier etwas posten zu können?

Bei mir hat's immer nur eine ganze weile gebraucht, wurde dann aber richtig übernommen.

[Ben04]

Jester schrieb:

sag bitte einfach bescheid wenn du kein feedback haben möchtest. ich lasse es dann. ansonsten verkneif dir bitte deinen imho ziemlich unverschämten ton.

Sorry war nicht so gemeint. Ich halte das Thema dennoch für zu komplex um es in Nebensätzen zu nennen.

Ich hab's mal weitergeleitet.

Danke

[Ben04]

Hab jetzt noch ein Beispiel zur Akkumulation hinzugefügt.

Ich wollte noch einen parallelen Quicksort hinzufügen, bei dem ich bei der naiven Implementierung anfangen würde. Das Problem ist, dass ich keine Erfahrungen habe, wie das besser geht, darum wollte ich mal fragen, ob das hier als "optimale" Version durchgeht? (Jester du hast das Beispiel doch schon mal angeführt. Kennst du dich da aus?)

template<class T>
class IsSmaller{
public:
	IsSmaller(const T&v):v(v){}
	bool operator()(const T&t)const{
		return t < v;	
	}
private:
	const T&v;
};

template<class T>
class IsEqual{
public:
	IsEqual(const T&v):v(v){}
	bool operator()(const T&t)const{
		return t == v;	
	}
private:
	const T&v;
};

template<class Iter>
void multithread_swap(Iter begin1, Iter end, Iter begin2, unsigned thread_count){
	int count = end - begin1;
	#pragma omp parallel for num_threads(thread_count)
	for(int i=0; i<count; ++i)
		std::swap(*(begin1+i), *(begin2+i));
}

template<class Iter, class Pred>
Iter multithread_partition(Iter begin, Iter end, Pred p, unsigned thread_count){
	assert(thread_count > 0);	
	if(thread_count == 1)
		return std::partition(begin, end, p);
	else{
		Iter mid = begin+(end-begin)/2;

		Iter cut1, cut2;

		#pragma omp parallel sections num_threads(2)
		{
			cut1 = multithread_partition(begin, mid, p, (thread_count)/2);
			#pragma omp section
			cut2 = multithread_partition(mid, end, p, (thread_count+1)/2);
		}	

		unsigned count = std::min(cut2-mid, mid-cut1);

		multithread_swap(cut1, cut1+count, cut2-count, thread_count);

		return cut1 + count;
	}
}

template<class T>
T median(const T&a, const T&b, const T&c){
	if(a <= b && b <= c)
		return b;
	else if(b <= a && a <= c)
		return a;
	else
		return c;
}

template<class Iter>
void multithread_quick_sort(Iter begin, Iter end, unsigned thread_count = omp_get_max_threads()){
	if(thread_count == 1)
		std::sort(begin, end); // Ich gehe davon aus, dass das ein sequentilles QuickSort ist
	else{
		if(end - begin <= 1)
			return;
		typedef typename std::iterator_traits<Iter>::value_type Value;

		Iter mid = begin+(end-begin)/2, last = end-1;
		Value pivot = median(*begin, *mid, *last);

		Iter 
			lower_bound = multithread_partition(begin, end, IsSmaller<Value>(pivot), thread_count),
			upper_bound = multithread_partition(lower_bound, end, IsEqual<Value>(pivot), thread_count);

		#pragma omp parallel sections num_threads(2) 
		{
			multithread_quick_sort(begin, lower_bound, thread_count/2);
			#pragma omp section
			multithread_quick_sort(upper_bound, end, (thread_count+1)/2);
		}
	}
}

[Ben04]

Seit längerer Zeit hab ich mal wieder ein wenig an diesem Artikel geschrieben. Was sich noch ändert:

• Eventuell ein paar Umformulierungen und/oder Verbesserung von Sprachschnitzern.
• Bei for Schleifen Teil werde ich noch auf die Neuerungen von OpenMP 3.0 eingehen. Diese bekommen allerdings nur einen Nebensatz da es noch einige Zeit dauern sollte ehe sie ausreichend verbreitet sind um produktiv eingesetzt zu werden.
• Bei den Beispielen kommt noch ein Beispiel für das Parallelisieren eines rekursiven Algorithmus rein. Wahrscheinlich das Quicksort Beispiel.

Der Rest wird gleich bleiben wenn es keine Verbesserungsvorschläge/Einwände gibt.

[Ben04]

Inhaltlich ist der Artikel nun fertig. Ich werde nichts mehr hinzufügen. Bitte meldet Teile die schwer zu verstehen, missverständlich oder gar falsch sind. Wenn niemand sich innerhalb der nächsten Woche meldet dann polier ich die Formatierung noch ein bisschen und setze den Artikel auf [R].

Der parallele Quicksort ist raus geflogen, da eine effiziente Implementierung doch trickreicher ist als ich dachte und den Rahmen sprengen würde. Zum Beispiel ist es bei der parallelen Version so wie ich sie oben vorgestellt habe wesentlich schlimmer wenn man das Pivot schlecht trift als bei der sequentiellen. Einmal schlecht getroffen und ein ganzer Prozessor ist quasi idle.

[Korbinian]

Hey, bin ausnahmsweise mal wieder hier aufgeschlagen weil ich das thema spontan interessiert hat - super Artikel, angenehm zu lesen und knackig erklaert! Wieder was gelernt

[Ben04]

Es sind noch 7 Tage bis zur nächsten Veröffentlichung und es hat sich noch kein Korrektor gemeldet. Ich möchte nicht stressen, nur weise ich lieber drauf hin falls es keiner gemerkt haben sollte, dass der Artikel auf [R] ist.

[darthdespotism]

~**Ich bin nun wirklich kein Rechtschreib-Genie daher hab' ich auch ein paar Stellen nicht weiter beachtet (v.a. Getrennt/Zusammenschreiben von Wörtern). Ansonsten ein toller Artikel

btw: Geht es nur mir so, dass die Smileys am Ende des Postings eingefügt werden und nicht an der Curserposition?**~

Inhaltsverzeichnis

• 1 OpenMP - Was ist das?

• 2 Compiler Unterstützung

• 3 Syntax

• 3.1 Parallel

• 3.2 Thread ID, single und master

• 3.3 Sections

• 3.4 Barrier

• 3.5 For

• 3.4.1 Verschachtelte Schleifen mit und ohne Collapse

• 3.4.2 Partielle sequentielle Abarbeitung mit ordered

• 3.6 Kurzschreibweise

• 4 Gemeinsamer Speicher

• 4.1 Ciritcal

• 4.2 Flush

• 4.3 Reduction

• 4.4 Locks/Mutex

• 5 Zeitmessung

• 6 Suchen in einem ungeordneten Array

• 7 Akkumulation

• 7.1 STL

• 8 Gibt es noch mehr?

• 9 Quellen

1 OpenMP - Was ist das?

OpenMP ist eine standardisierte Schnittstelle für Fortran, C und C++, welche die Entwickelung von Programmen, welche mehrere Prozessor-Kerne nutzen, vereinfachen soll. Ich werde in diesem Artikel aber nur auf die C++ Schnittstelle eingehen.

2 Compiler Unterstützung

Es gibt mehrere Versionen von OpenMP. Dieser Artikel beschäftigt sich hauptsächlich mit OpenMP 2.5. Ich werde allerdings auch auf ein paar Neuerungen in der Version 3.0 eingehen. OpenMP wird von einer Reihe von Compilern unterstützt, jedoch meistens nicht in den Basisausstattungen dieser. Mit dem GCC gibt es aber mindestens ein Compiler der jedem Leser zugänglich ist.

Folgende Tabelle liefert einen kleinen Überblick.

• GCC >= 4.2 (mit dem -fopenmp Parameter)
• MinGW Technology-Preview-Version (mit dem -fopenmp Parameter)
• VC 2005 Professional und Team System (mit dem /openmp Parameter, Standard und Express Versionen gehen nicht)

Auf der OpenMP-Webseite findet man eine detailliertere Tabelle.

3 Syntax

OpenMP besteht zu einem Großteil aus Compiler-Erweiterungen. Es gibt eine Reihe von Pragma-Anweisungen welche die OpenMP-Funktionalität bereitstellen. Jede von ihnen hat folgenden Syntax:

#pragma omp Anweisung
C++ Code auf den sich die Anweisung bezieht

Neben diesen pragma-Konstrukten gibt es noch eine Reihe von normalen Funktionen. Diese werden im Header omp.h definiert. Alle haben den Präfix omp_.

Jeder Compiler der sich an den C++-Standard hält muss unbekannte Pragma-Anweisungen ignorieren und höchstens mit einer Warnung versehenen. Fehler sich nicht standard-konform. Durch das Benutzen eines Dummy-Headers, ist es deswegen in vielen Fällen möglich eine OpenMP-Anwendung sequentiell auf nicht OpenMP-fähigen Compilern zu übersetzen. Dies ist auch empfehlenswert wenn der Compiler zwar OpenMP unterstützt, allerdings der Zielprozessor keine echte Parallelisierung unterstützt.

In diesem Artikel werde ich versuchen die Anweisungen welche kein sequentielles Fallback-Konstrukt haben zu ignorieren. Sie bieten auch eher selten Funktionalität an, welche man nicht auch anders erhalten kann. Daher ist dies an sich kein großer Verlust.

3.1 Parallel

Die grundlegend Anweisung, ohne die bei OpenMP nichts auskommt, ist die parallel-Anweisung. Sie erzeugt eine Reihe von Threads und führt den nachfolgenden C++ Block mit jedem Thread aus. Dies bedeutet dass der gleiche Code mehrmals gleichzeitig ausgeführt werden kann und zwar von jedem Thread einmal. Wie viele Threads erstellt werden, ist der Implementierung von OpenMP freigestellt. Meistens werden aber so viele Threads erzeugt wie es im Rechner Kerne gibt.

Es ist also in der Regel nicht möglich vorher zu sagen wie oft ein parallel-Block ausgeführt wird.

int main(){
	printf("Hier ist noch alles einfach und sequentiell.\n");
	#pragma omp parallel
	{
		printf("Hallo ");
		printf("Welt!\n");
	}
	printf("Alle Threads haben Hallo gesagt.");
}

Im parallel-Block ist ganz normaler C++-Code erlaubt. Man kann C++-Objekte anlegen und auch Funktionen aufrufen. Sogar Ausnahmen sind erlaubt, solange sie nicht aus dem parallel-Block rausfliegen.

Solange sich die Threads nicht in die Quere kommen, ist das auch schon alles was man machen muss um sämtliche Kerne seines Rechners zu benutzen. Allerdings ist dies leider nur sehr selten der Fall~ein ist zu viel~. Bereits in diesem einfachen Fall ist es möglich, dass Threads sich streiten. Wenn zwei Threads gleichzeitig etwas ausgeben, ist es nämlich nicht klar welcher zuerst dran kommt. Alles von Funktionieren-wie-erwartet bis hin zu Abstürzen ist eine möglich Konsequenz.

Was das Problem noch komplizierter macht, ist dass das Verhalten nicht vollständig deterministisch ist. Wie viele Threads erzeugt werden kann von Lauf zu Lauf unterschiedlich sein und welcher Thread wann dran kommt kann man noch schlechter vorhersagen. Dies führt dazu, dass manche Bugs nur mit einer gewissen Wahrscheinlichkeit eintreten. Das Debuggen solcher Fehler ist sehr schwer.

An sich müsste ich hier sicherstellen, dass zu jedem Zeitpunkt immer höchstens nur ein Thread etwas auf die Konsole schreibt. Bei printf kann man aber normalerweise davon ausgehen, dass legendlich die Reihenfolge der Ausgaben etwas durcheinander geraten kann. Abstürze gehören eher in den Bereich der Theorie. printf leidet normalerweise weniger unter diesem Problem als std::cout, deswegen werde ich in diesem Artikel ausnahmsweise printf bevorzugen. Im Folgenden werde ich davon ausgehen, dass Ausgaben mit printf einfach funktionieren. Alles andere würde die Beispiele zu kompliziert machen.

Da man nicht wirklich weiß wie oft ein parallel-Block ausgeführt wird, ist eine Parallelisierung recht umständlich ohne weitere OpenMP-Konstrukte. Er eignet sich aber in Regel vorzüglich um threadspezifische Variablen zu initialisieren.

Mit Hilfe der num_threads(n)-Anweisung kann man festlegen wie viele Threads erstellt werden. n muss eine Ganzzahl sein und braucht nicht zum Zeitpunkt der Übersetzung bekannt zu sein.

int main(){
	printf("Wie viele Threads sollen es sein?\n");
	int n;
	scanf("%d", &n);
	#pragma omp parallel num_threads(n)
	{
		printf("Hallo ");
		printf("Welt!\n");
	}
	printf("Alle Threads haben Hallo gesagt.");
}

Wenn man für n eine große Zahl eingibt, kann man mit großer Wahrscheinlichkeit beobachten wie sich die Threads in die Quere kommen. Es folgt eine Beispielausgabe, so wie ich sie beobachten konnte.

Hallo Hallo Welt!
Hallo Welt!
Welt!
Hallo Welt!

Während der eine Thread noch am Schreiben ist, funkt ein anderer ihm einfach dazwischen. Versuchen Sie nach zu vollziehen wie diese Ausgabe zustande kam. Es ist sehr wichtig dieses Problem zu verstehen, denn sonst werden Sie mit OpenMP oder auch Multithread-Programmierung im Allgemeinen nur Ärger haben.

omp_get_max_threads gibt an wie viele Threads im nächsten parallel-Block erzeugt werden.

Bei den meisten OpenMP-Implementierungen werden die Threads nur einmal erstellt. Dies gilt auch wenn es mehrere parallel-Blocks gibt. Zwischen den Blocks gehen die zusätzlichen Threads in einen Wartezustand. Sie werden erst am Ende des Programms zerstört. Die in Regel recht hohen Threaderstellungskosten muss man also nur genau einmal zahlen.

parallel-Blöcke dürfen verschachtelt werden. In der Standardeinstellung erzeugt der verschachtelte Block aber keine neuen Threads.

3.2 Thread ID, single und master

Jeder Thread besitzt eine eindeutige Identifikationsnummer. Es handelt sich dabei um fortlaufende Nummern und nicht etwa um gecastete Zeiger oder ähnliches. Man kann die ID also als Index in einem Array benutzen. Den Thread mit der ID 0 nennt man auch noch Master.

Mit Hilfe von omp_get_thread_num kann man die ID eines Threads bestimmen und durch einen Aufruf von omp_get_num_threads kann man die Gesamtanzahl der Threads, welche vom parallel-Block erzeugt wurden, heraus finden.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());	
		if(omp_get_thread_num() == 0)
			printf("Es gibt %d Threads.\n", omp_get_num_threads());		
	}

}

Wie man an diesem Beispiel sieht, ist es möglich auch in einem parallel-Block, Instruktionen nur einmal aus zu führen. Allerdings gibt es dafür bessere OpenMP-Konstrukte.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());	

		#pragma omp master
		{
			printf("Es gibt %d Threads.\n", omp_get_num_threads());	
		}

		#pragma omp single
		{
			printf("Thread %d war hier.\n", omp_get_thread_num());	
		}

		#pragma omp single
		{
			printf("Danach war Thread %d da.\n", omp_get_thread_num());	
		}
	}
}

Ein master-Block wird nur vom Master-Thread ausgeführt. Der Code ist äquivalent zur if-Abfrage aus dem vorherigen Beispiel. Ein single-Block wird genau einmal ausgeführt, allerdings ist es egal von welchem Thread.

Sämtliche Threads warten am Ende eines single-Blocks bis auch der letzte es bis dahin geschafft hat. Dadurch ist sichergestellt, dass der zweite single Block nicht ausgeführt wird ehe der erste vollständig abgearbeitet wurde. Durch eine nowait Anweisung kann man das Warten unterdrücken.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			printf("Manchmal bin ich als erstes da.\n");	
		}

		#pragma omp single nowait
		{
			printf("Hin und wieder gewinne aber auch ich das Rennen.\n");	
		}
	}
}

Wenn man nowait verwendet dann wird nur noch garantiert, dass die Ausführung des ersten Blocks vor der des zweiten beginnt. Es ist aber gut möglich, dass der zweite vor dem ersten beendet wird. Es ist auch möglich, dass die Ausführung des ersten Threads noch vor der ersten Anweisung im Block unterbrochen wird und es dadurch erscheint, als wenn der zweite Block vor dem ersten gestartet werden würde. In der Praxis kann man also nicht darauf bauen, dass der erste vor dem zweiten gestartet wird, auch wenn dies mit einer hohen Wahrscheinlichkeit so sein wird.

Am Ende eines Master-Block wird nicht gewartet.

Bei verschachtelten parallel-Blöcken werden die IDs durch den innersten Block neu verteilt. In der Standardeinstellung erhalten alle Threads die ID 0 und die Anzahl der Threads wird als 1 gemeldet. Die Idee dahinter ist, dass sich alles auf den innersten Block bezieht. Dieser hat aber keinen Thread erzeugt und folglich gibt es nur den der den Block angestoßen hat. Eine weitere Konsequenz ist, dass single- und master-Blöcke mehrfach ausgeführt werden. Ein verschachtelter parallel-Block wird demnach also in der Standardeinstellung sequentiell abgearbeitet. Mit num_threads(n) kann das erzeugen weiterer Threads erzwingen.

Die Problematik wird durch folgendes Beispiel gut illustriert.

int main(){
	#pragma omp parallel
	{
		#pragma omp parallel
		{
			printf("omp_get_thread_num() = %d,"
				" omp_get_num_threads() = %d\n", 
				omp_get_thread_num(), 
				omp_get_num_threads());	
		}
	}
}

3.3 Sections

Alternativ zu den single-Blöcken kann man auch sections verwenden.

int main(){
	#pragma omp parallel
	{
		#pragma omp sections
		{
			{
				printf("Thread %d war hier.\n", omp_get_thread_num());			
			}
			#pragma omp section
			{
				printf("Aber hier war Thread %d.\n", omp_get_thread_num());		
			}
			#pragma omp section
			{
				printf("und eine weitere section");		
			}
		}
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

Jede der drei Sektionen wird genau einmal ausgeführt. Am Ende des sections-Block warten die Threads ehe sie fortfahren. Auch hier kann man das Warten durch nowait unterdrücken.

Man könnte obiges Beispiel auch mit single ausdrücken.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			printf("Thread %d war hier.\n", omp_get_thread_num());			
		}
		#pragma omp single nowait
		{
			printf("Aber hier war Thread %d.\n", omp_get_thread_num());		
		}
		#pragma omp single nowait
		{
			printf("und eine weitere section");		
		}
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

Bei der single-Variante ist der erste Block mit hoher Wahrscheinlichkeit der erste der gestartet wird. Bei der sections-Variante sind alle Blöcke gleich wahrscheinlich.

3.4 Barrier

Wenn an einer beliebigen Stelle gewartet werden soll so kann man eine barrier-Anweisung benutzen.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());
		#pragma omp barrier
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

In diesem Beispiel ist sicher gestellt, dass alle Threads "Hallo" sagen ehe sich der erste verabschiedet.

3.5 For

Mit Hilfe von OpenMP kann man auch einfache Schleifen parallelisieren.

int main(){
	#pragma omp parallel
	{
		#pragma omp for
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);
	}
}

Anders als bei einer normalen Schleife, werden die Durchläufe hier nicht unbedingt nacheinander abgearbeitet, sondern möglicherweise gleichzeitig und auch in zufälliger Reihenfolge. Dies führt dazu, dass die Ausgabe des obigen Codes nicht immer die gleiche ist. Die Reihenfolge der Zeilen kann beliebig permutiert sein.

Um dies zu bewerkstelligen werden die Durchläufe möglichst gleichmäßig zwischen den Threads verteilt. Zum Beispiel bei zwei Threads, könnte der erste die Durchläufe mit 0<=i<=4 abarbeiten und der andere die restlichen mit 5<=i<=9. Die genaue Verteilung ist dem Compiler frei gestellt. Sie hängt jedoch nur von der Anzahl der Threads und der Anzahl der Schleifendurchläufen ab, das heißt wenn aus irgendeinem Grund ein Thread massive Verspätung hat, dann wird nicht kurzerhand die Arbeit umgeschichtet.

int main(){
	#pragma omp parallel
	{
		// Sorgt für massive Verspätung eines Threads
		#pragma omp single nowait
		{
			for(int i=0; i<INT_MAX; ++i)
				for(int j=0; j<INT_MAX; ++j)
					{}
		}

		// Tja hier warten dann die anderen...
		#pragma omp for
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);
	}
}

Obwohl der eine Thread sehr lange mit dem single-Block beschäftigt ist, kommt keiner der wartenden Threads auf die Idee bei der for-Schleife Überstunden zu leisten. Dies ist nur mit einer dynamischen Arbeitsverteilung möglich.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			for(int i=0; i<INT_MAX; ++i)
				for(int j=0; j<INT_MAX; ++j)
					{}
		}

		#pragma omp for schedule(dynamic)
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);

		// Hier wird erst gewartet wenn die Schleife vollständig abgearbeitet wurde.
	}
}

Wenn ein Thread fertig ist, schaut er ob es noch Arbeit gibt und übernimmt diese dann auch gleich. Ein ähnliches Problem ergibt sich wenn die Durchläufe unterschiedlich lang dauern.

Am Ende eines for-Blocks wird gewartet. Dies kann man wie üblich mit nowait unterdrücken.

Da schedule(dynamic) einen größere Verwaltungskosten als die normale Version hat, kann es sein, dass die normale Version in manchen Fällen schneller ist, obwohl nicht alle Threads immer beschäftigt sind. Die normale Version kommt nämlich ganz ohne Kontakt zwischen den einzelnen Threads zurecht und ist daher sehr schnell. Ich werde weiter unten auf konkrete Implementierungsmöglichkeiten eingehen.

Damit eine Schleife überhaupt parallelisiert werden kann muss sie eine Reihe von sehr strikten Kriterien erfüllen.

• Die Schleifenvariable muss eine Ganzzahl sein. In OpenMP 3.0 sind auch Random-Access-Iteratoren und Zeiger erlaubt.
• Die Abbruchbedingung muss ein Ordnungsvergleich (<,>, <= oder sein und auf der linken Seite muss sich die Schleifenvariable und nur die befinden. Die rechte Seite darf nicht während der Ausführung der Schleife verändert werden.
• Die Inkrementierung muss in immer gleich großen Schritten stattfinden. Dies bedeutet, dass man nur --, ++, += oder -= verwenden darf um die Schleifenvariable zu verändern. Die Größe der Schritte darf nicht während der Ausführung der Schleife verändert werden.
• break und return sind verboten. continue ist erlaubt. goto darf nur verwendet werden wenn man nicht aus der Schleife raus springt.

Folgende Schleifen lassen sich nicht mit OpenMP parallelisieren.

int n = foo();

// Auf der einen Seite der Bedingung muss ein einfaches i stehen.
for(int i=0; i*i<n; ++i)
	printf("%d", i);

// Die Inkrementierungsschritte sind nicht gleichgroß
for(int i=0; i<n; i*=3)
	printf("%d", i);

// Dies ist erlaubt, allerdings ist unspezifiziert wie oft bar und foobar aufgerufen werden.
for(int i=-10; i<n*n*bar(); i+=foobar(n))
	printf("%d", i);

// Man darf nicht aus der Schleife raus springen.
int a[10] = { ... }
for(int j=0; j<10; ++j)
	if(a[j] == n)
		return i;

// n darf nicht in der Scheife verändert werden.
for(int i=0; i<n*10; i+=n)
	++n;

3.4.1 Verschachtelte Schleifen mit und ohne Collapse

OpenMP-for-Schleifen dürfen nicht verschachtelt werden. Man kann allerdings mit einem zusätzlichen parallel-Block Abhilfe schaffen. Dies ist jedoch selten sinnvoll da die Schleife zwar eine OpenMP-Schleife ist, aber dennoch sequentiell abgearbeitet wird.

#pragma omp parallel
{
	#pragma omp for
	for(int i=0; i<10; ++i)
	{
		#pragma omp parallel
		{
			#pragma omp for
			for(int j=0; j<10; ++j)
			{
				printf("%d*%d = %d\n", i, j, i*j);
			}
		}
	}
}

Es wird aber nur die äußere Schleife parallelisiert und die innere sorgt nur für Probleme und Kosten. Besser ist:

#pragma omp parallel
{
	#pragma omp for
	for(int i=0; i<10; ++i)
	{
		// Dieses for ist ein ganz normales C/C++ for
		for(int j=0; j<10; ++j)
		{
			printf("%d*%d = %d\n", i, j, i*j);
		}

	}
}

In OpenMP 3.0 gibt es eine spezielle collapse-Schleife um über ein mehrdimensionales Array zu iterieren. Dies entspricht einem der häufigsten Einatzgebiet von mehrfach verschachtelten Schleifen. Die Schleifen müssen allerdings sehr einfach gehalten sein, damit eine Parallelisierung überhaupt erfolgen kann. Alle Schleifenvariablen müssen völlig unabhängig von einander sein. Es darf nur die innerste Schleife beliebigen Code enthalten. Alle anderen dürfen nur verschachtelte Schleifen enthalten. Das Argument der collapse-Anweisung gibt an wie tief die Schleifen verschachtelt sind.

Dies sind sehr starke Einschränkungen und deswegen deckt das folgende Beispiel an sich auch quasi jede erlaubte Variante einer collapse-Schleife ab.

int main(){
	#pragma omp parallel
	{
		int foo[10][3][7];
		#pragma omp for collapse(3)
		for(int i=0; i<10; ++i)
			for(int j=0; j<3; ++j)
				for(int k=0; k<7; ++k){
					// Nur hier darf beliebiger Code stehen!
					foo[j][k] = 0;
				}

}

3.4.2 Partielle sequentielle Abarbeitung mit ordered

Man kann in eine OpenMP-[i]for* Schleife einen ordered-Block packen. Dieser sorgt dafür, dass ein Teil der Schleife sequentiell abgearbeitet wird.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<=10; ++i){
			double v = i;
			for(int j=1; j<10; ++j)
				v *= i;
			#pragma omp ordered
			{
				printf("%d^10 = %.0f\n", i, v);
			}
		}
	}
}

Die Berechnung der Werte kann in beliebiger Reihenfolge und parallel erfolgen. Die Ausgabe ist jedoch sequentiell, das heißt die Ausgabe ist immer die gleiche.

0^10 = 0
1^10 = 1
2^10 = 1024
3^10 = 59049
4^10 = 1048576
5^10 = 9765625
6^10 = 60466176
7^10 = 282475249
8^10 = 1073741824
9^10 = 3486784401
10^10 = 10000000000

Ein ordered-Block darf nur in einer speziellen ordered-for-Schleife auftreten. Jeder Durchlauf darf nur maximal einem ordered-Block begegnen. Folgendes Beispiel ist demnach falsch.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<10; ++i){
			#pragma omp ordered
			printf("%d^10 =", i);

			double v = i;
			for(int j=1; j<10; ++j)
				v *= i;

			// Hier wird etwas schief gehen da 
			// dies der zweite Block ist.
			#pragma omp ordered
			printf("%.0f\n", v);
		}
	}
}

Dies bedeutet jedoch nicht, dass man nur einen ordered-Block pro Schleife haben darf.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<10; ++i){
			if(i%3 == 0){
				#pragma omp ordered
				printf("%d ist ein Vielfaches von 3\n", i);
			}else{
				#pragma omp ordered
				printf("%d ist kein ein Vielfaches von 3\n", i);
			}
		}
	}
}

In dem Beispiel ist sichergestellt, dass kein Durchlauf zwei verschiedene Blöcke zu sehen bekommt.

3.6 Kurzschreibweise

Wenn ein parallel-Block nur eine for-Schleife enthält, dann gibt es eine Kurzschreibweise.

#pragma omp parallel for
for(int i=0; i<10; ++i)
{
	for(int j=0; j<10; ++j)
	{
		printf("%d*%d = %d\n", i, j, i*j);
	}

}

Analog gibt es auch ein parallel sections-Block.

Folgendes Beispiel berechnet alle Primzahlen kleiner als eine gegebene Zahl.

bool is_prime(unsigned n){
	for(unsigned i=2; i<n; ++i)
		if(n%i == 0)
			return false;
	return true;
}

int main(){
	unsigned n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%u", &n);

	// schedule(dynamic) ist notwendig da is_prime keine
	// konstante Laufzeit besitzt. 
	#pragma omp parallel for schedule(dynamic)
	for(int i=2; i<=n; ++i)
		if(is_prime(i))
			printf("%u ist prim\n", i);
}

4 Gemeinsamer Speicher

Die Zugriffsregeln für Variablen sind die die man erwarten würde. Variablen leben so lange wie der Block in dem sie definiert wurden. Sie können verdeckt werden und man kann auf Variablen einer höheren Ebene zugreifen.

int a;

int main(){
	int b;
	#pragma omp parallel
	{
		int c;
		// a und b werden geteilt, jeder Thread hat sein eigenens c
	}
}

Wenn ein Thread schreibend auf a oder b zugreift dann ergeben sich sämtliche vom Multithreading her bekannte Probleme. OpenMP bietet aber leider nur wenig Revolutionäres um dieses alt bekannte Problem zu lösen.

4.1 Ciritcal

Ein critical-Block ist ein Block in dem sich immer nur genau ein Thread befindet. Wenn bereits ein Thread im Block ist, wartet jeder Thread der hinein will, bis der Block wieder frei wird.

Betrachten Sie folgendes Beispiel:

int main(){
	#pragma omp parallel for
	for(int i=0; i<10; ++i){
		int n;
		#pragma omp critical
		{
			printf("Geben sie eine Zahl ein: ");
			scanf("%d", &n);
		}
		printf("%d * %d = %d\n", i, n, i*n);
	}
}

Hier wird ein critical-Block verwendet um sicher zu stellen, dass der Benutzer nie aufgefordert wird zwei oder mehr Zahlen gleichzeitig einzugeben. Nimmt man das critical weg dann kann es zu folgender Ausgabe kommen.

Geben sie eine Zahl ein: Geben sie eine Zahl ein:

Jeder critical-Block besitzt einen Namen. Wenn zwei Blöcke denselben Namen besitzen, kann nur genau ein Thread in einem der beiden sein. Zwei Blöcke mit verschiedenen Namen beeinflussen sich nicht. Das heißt es ist durchaus möglich, dass 2 Threads in critical-Blöcken sind sofern sie verschiedene Namen haben.

Wenn bei einem critical-Block kein Name angegeben wird, so wie im ersten Beispiel, besteht der Name aus dem leeren Wort. Alle leeren Wörter werden als gleich angesehen. Im Klartext bedeutet dies, dass keine zwei Threads gleichzeitig in irgendeinem namenlosen Block sein können.

Leider sind die Namen global und liegen in keinem Namensraum. Dadurch sind Namenskonflikte möglich. Es ist allerdings sehr unwahrscheinlich, dass dies andere Auswirkungen hat als einfach nur ein langsameres Programm, da lediglich ein Thread warten muss obwohl dies nicht nötig gewesen wäre. Es ist sichergestellt, dass sich normalen Symbolen wie Funktionsnamen nicht in die Quere kommen.

int main(){
	FILE*file = fopen("out.txt", "w");
	#pragma omp parallel for
	for(int i=0; i<10; ++i){
		int n;

		#pragma omp critical(console)
		{
			printf("Geben sie eine Zahl ein : ");
			scanf("%d", &n);
		}

		#pragma omp critical(console)
		{
			printf("%d * %d = %d\n", n, n, n*n);
		}

		#pragma omp critical(file)
		{
			fprintf(file, "%d * %d = %d\n", n, n, n*n);
			fprintf(file, "%d + %d = %d\n", n, n, n+n);
		}

		#pragma omp critical(console)
		{
			printf("%d + %d = %d\n", n, n, n+n);
		}
	}
	fclose(file);
}

In diesem Beispiel ist sicher gestellt, dass nichts auf die Konsole geschrieben wird, falls der Benutzer aufgefordert wurde etwas einzugeben. In der Datei müssen beide zueinander passenden Rechnungen hinter einander ausgegeben werden. Auf der Konsole ist dies nicht garantiert.

4.2 Flush

Viele Rechner besitzen Register und mehrere Level von Speicher-Caches, um den Zugriff auf Variablen zu beschleunigen. Anstatt direkt mit der Variable im Hauptspeicher zu arbeiten wird diese zuerst in einen schnelleren Cache-Speicher kopiert. Anschließend wird mit der Kopie gerechnet und nach einiger Zeit wird die Variable in den Hauptspeicher zurück kopiert. Dies führt dazu, dass es mehrere Versionen von der gleichen Variable gibt welche nicht immer gleich sind.

In einem sequentiellen Programm stellt dies kein Problem dar, da zu jedem Zeitpunkt bekannt ist wo sich die rezenteste ~aktuelleste ist meiner meinung nach einfacher zu verstehen - das Wort ist mir noch nie begegnet~ Version der Variable befindet. Jedoch ist dies bei mehreren Threads nicht mehr so einfach~Yoda Speek~. Der eine weiß nicht wann der andere fertig ist und wann er die Werte synchronisieren muss.

Um dieses Problem zu lösen bietet OpenMP die flush-Anweisung an. Sie synchronisiert die lokale Version des Threads mit der des Hauptspeichers. Dies bedeutet, dass jeder Thread nachdem er fertig ist, eine geteilte Variable zu verändern, diese auch flushen muss. Ein Thread der lesend auf eine Variable zugreift will, muss die Variable auch flushen damit er sicher ist die neuste Version zu Gesicht zu bekommen. Wenn ein Thread auch mit einer veralteten Version arbeiten kann, braucht er die Variable nicht zu flushen. Wenn sicher ist, dass eine Variable nicht in einem parallel-Block verändert wird, dann kann man sich das flush auch sparen da es in alle Versionen gleich sind, ganz egal wo sie sich befinden. ~Besser vll: da alle Threads die gleiche Version sehen, ...~

Folgendes Programm berechnet das Faktoriell einer Zahl.

int main(){
	unsigned n, fact = 1;
	printf("Geben sie eine Zahl ein: ");
	scanf("%u", &n);

	#pragma omp parallel
	{
		unsigned local_fact = 1;

		// Auf n kann man ohne flush zugreifen,
		// da es nicht verändert wird.
		#pragma omp for
		for(int i=1; i<=n; ++i){
			// local_fact braucht nicht geflusht 
			// zu werden, da es nicht geteilt wird.
			local_fact *= i;
		}

		#pragma omp critical
		{
			// Zuerst holt der Thread sich die neuste Version
			// aus dem Hauptspeicher.
			#pragma omp flush(fact)

			// Danach wird die Variable verändert.
			fact *= local_fact;

			// Anschließend wird die Variable wieder in den
			// Hauptspeicher geschrieben, damit andere Threads
			// sich die neuste Version von dort laden können.
			#pragma omp flush(fact)
		}	
	}

	printf("%u! = %u\n", n, fact);
}

Ein fehlendes flush bedeutet nicht, dass die Variablen nicht synchronisiert werden. Es steht der Implementierung frei dies zu tun wann immer sie will. Dies führt leider oft zu nur sehr schwer reproduzierbaren Bugs. Zum Beispiel ist es sehr wahrscheinlich, dass das Betriebssystem die Rechnerzeit anders auf die Threads verteilt, wenn die Hardware unter voller Last läuft, als wenn es nur mit kleinen Texteingaben konfrontiert wird. Es ist daher sinnvoll im Zweifelsfall lieber ein flush zuviel zu schreiben.

4.3 Reduction

Mit reduction lassen sich Berechnung vom Typ s += f(1)+f(2)+...+f(n) parallelisieren. Jeder Thread rechnet ein paar Terme aus und summiert die auf welche er berechnet hat. Nachdem alle Threads fertig sind zählt jeder Thread sein Zwischenergebnis zur Variable s hinzu.

Diese Zwischenergebnis wird in einer Variable mit dem gleichen Namen wie s gespeichert. Diese Hilfsvariable wird im reduction-Block angelegt und mit dem neutralen Element (0 im Fall der Addition) initialisiert. Man sollte sich aber nicht auf diese***s*** Verhalten verlassen. Wenn ein Programm sequentiell übersetzt wird, das heißt ohne OpenMP Unterstützung, dann gibt es keine Zwischenergebnisse und deswegen hat s beim ersten Durchlauf auch nicht unbedingt den Wert des neutralen Elements.

Neben + lassen sich auch -, *, &, |, ^, && und || auf diese Weise parallelisieren. Es muss bei den Operatoren jedoch um ihre Buildin-Version handeln.

• +, -, | und ^ besitzen das neutrale Element 0.
• * besitzt das neutrale Element 1.
• & besitzt das neutrale Element ~0. (Bei ~0 sind alle Bits gesetzt sind.)
• && besitzt das neutrale Element true.
• || besitzt das neutrale Element false.

Bei Gleitkommazahlen hängt die Genauigkeit von der Auswertungsreihenfolge ab. Die Reihenfolge ist nicht spezifisiert, also auch nicht die Genauigkeit des Ergebnisses. Es ist nicht einmal garantiert, dass es immer das gleiche Ergebnis sein wird. Gleitkommazahlen in Kombination mit einer Reduktion haben einen gewissen nicht deterministischen Teil im Ergebnis.

In der Fort***r***an Version von OpenMP gibt es zusätzlich noch Min- und Max-Reduktionen. Die wurden aus mir unerklärlichen Gründen in der C/C++ Version vergessen.

Folgendes Beispiel bestimmt ob eine Zahl prim ist.

int main(){
	int n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%d", &n);

	// Der Initialwert von is_prime ist wichtig da dieser auch
	// in die Berechnung einbezogen wird.
	bool is_prime = true;

	#pragma omp parallel reduction(&&: is_prime)
	{
		// jeder bekommt sein eigenes is_prime,
		// welches mit true, dem neutralen Element von &&,
		// initialisiert wurde.

		#pragma omp for
		for(int i=2; i<n; ++i)
			if(n%i == 0)
				is_prime = false;

		// Alle is_prime (auch das aus der main Funktion) werden 
		// mittels && verknüpft. Im Klartext bedeutet dies, dass 
		// falls eines false ist, ist das Resultat false.
	}

	if(is_prime)
		printf("%d ist prim\n", n);
	else
		printf("%d ist nicht prim\n", n);

}

Am Ende des parallel Blocks wird folgende Rechnung ausgeführt:

is_prime = is_prime && thread[0].is_prime && thread[1].is_prime && ...  && thread[thread_num-1].is_prime

Bei der thread[j]-Notation handelt es sich legendlich um Pseudocode.

Auch hier gibt es eine Kurzschreibweise. Folgendes Beispiel berechnet das Faktoriell einer Zahl.

int main(){
	unsigned n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%u", &n);

	unsigned fact = 1;

	#pragma omp parallel for reduction(*: fact)
	for(int i=2; i<n; ++i)
		fact *= i;

	printf("%u! = %u\n", n, fact);
}

reduction erlaubt es dem Compiler eine Reihe von Optimierungen durchzuführen, welche bei der äquivalenten Formulierungen mit critical und flush nur sehr schwer möglich sind.

4.4 Locks/Mutex

OpenMP bietet auch Mutexe an. Hier weicht OpenMP an sich nicht von allgemein bekanten Konzepten ab. Es gibt einen Lock-Type namens omp_lock_t. Jeder Mutex ist entweder offen oder geschlossen. Folgende Funktionen werden bereitgestellt um damit zu arbeiten:

• void omp_init_lock(omp_lock_t*)
  Initialisiert den Mutex. Der Mutex ist danach offen.
• void omp_destroy_lock(omp_lock_t*)
  Zerstört den Mutex.
• void omp_set_lock(omp_lock_t*)
  Versucht den Mutex zu schließen. Wurde der Mutex bereits von einem anderen Thread geschlossen so wird gewartet bis, dass dieser ihn wieder geöffnet hat.
• void omp_unset_lock(omp_lock_t*)
  Öffnet den Mutex wieder
• int omp_test_lock(omp_lock_t*)
  Versucht den Mutex zu schließen. Wurde der Mutex bereits von einem anderen Thread geschlossen so wird 0 zurück gegeben, ansonsten 1.

omp_lock_t ist ein nicht rekursives Mutex. Das heißt ein Thread darf einen Mutex nicht mehrmals schließen. Er würde in dem Fall nämlich auf sich selbst warten. Es gibt aber auch eine rekursive Variante und die heißt omp_nest_lock_t. Die Funktionen um damit umzugehen heißen ähnlich und funktionieren analog.

• void omp_init_nest_lock(omp_next_lock_t*)
• void omp_destroy_nest_lock(omp_next_lock_t*)
• void omp_set_nest_lock(omp_next_lock_t*)
• void omp_unset_nest_lock(omp_next_lock_t*)
• int omp_test_nest_lock(omp_next_lock_t*)

Da diese Schnittstelle alles anderes als modernes C++ ist, bietet es sich an diese Funktionen zu wrappen.

class Mutex
{
public:
	Mutex()				{omp_init_lock(&lock);}
	~Mutex()			{omp_destroy_lock(&lock);}
	void lock()			{omp_set_lock(&lock);}
	void unlock()			{omp_unset_lock(&lock);}
	bool try_to_lock()		{return omp_test_lock(&lock);}
private:
	Mutex(const Mutex&);
	Mutex&operator=(const Mutex&);
	mp_lock_t lock;
};

class ScopedLock
{
public:
	ScopedLock(Mutex&lock)		:lock(lock){lock.lock();}
	~ScopedLock()			{lock.unlock();}
private:
	ScopedLock(const ScopedLock&);
	ScopedLock&operator=(const ScopedLock&);
	Mutex&lock;
};

5 Zeitmessung

Neben den Funktionen welche direkt mit der Multithread-Programmierung bietet OpenMP noch zwei nützliche Funktionen zur Zeitmessung an.

• double omp_get_wtime();
  Gibt die momentane Zeit in Sekunden zurück. Es ist nicht definiert wann der Zeitpunkt 0 ist. Es ist legendlich bekannt, dass es nicht während dem Lauf eines Programms verändert wird und in der Vergangenheit liegt. Dies ist aber nicht weiter schlimm, wenn man die Funktion nur für Zeitmessungen einsetzen will. Die Genauigkeit dieser Funktion ist nicht spezifisiert und darf leicht von Thread zu Thread abweichen. Man kann die Genauigkeit aber mit omp_get_wtick in Erfahrung bringen.
• double omp_get_wtick();
  Gibt den kleinstmöglichen Abstand zwischen 2 verschiedenen von omp_get_wtime zurückgegebenen Zeitpunkten an. Der Rückgabewert von omp_get_wtick darf von Thread zu Thread variieren.

6 Suchen in einem ungeordneten Array

Das Suchen in einem unsortierten Array lässt sich sehr leicht mit OpenMP beschleunigen.

// Ohne OpenMP
int find_first(int arr[], int size, int value){
	for(int j=0; j<size; ++j)
		if(arr[j] == value)
			return j;
	return -1;
}

// Mit OpenMP
int find_any(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel for
	for(int j=0; j<size; ++j)
		if(arr[j] == value){
			i = j;
			#pragma omp flush(i)		
		}
	return i;
}

Das aktuelle i braucht nicht zuerst per flush-Anweisung geladen zu werden, da der Wert eh nur überschrieben werden würde.

Wie man bereits am Namen erkennt, ist bei der parallelen Version nicht definiert welchen Wert gefunden wird, wenn es mehrere gleiche Werte im Array gibt. Wenn dies jedoch nötig ist, kann man dieses Problem leicht mit ordered beheben.

int find_first(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel for ordered
	for(int j=size-1; j>=0; --j)
		if(arr[j] == value){
			#pragma omp ordered
			i = j;
			#pragma omp flush(i)		
		}
	return i;
}

Eines der Nachteile der parallelen Version ist, dass immer alle Elemente angeschaut werden. Dies ist aber nicht immer notwendig denn man kann die Schleife abbrechen sobald ein Element gefunden wurde. OpenMP bietet allerdings keine Möglichkeit dies zu tun. Man kann versuchen die Schleife nach zu bauen.

int find_any(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel
	{
		const int begin = size*omp_get_thread_num() / omp_get_num_threads();
		const int end = size*(omp_get_thread_num()+1) / omp_get_num_threads();

		for(int j=begin; j<end; ++j){
			if(arr[j] == value)
				i = j;
			#pragma omp flush(i)
			if(i != -1)
				break;
		}
	}
	return i;
}

Dieser Code dürfte recht nahe an dem sein, was der Compiler aus einer OpenMP-Schleife mit statischer Verteilung macht. In der Regel sollte man aber die OpenMP-for-Schleife bevorzugen. Diese ist weit weniger fehleranfällig. So kann es beim Beispiel oben zum Beispiel zu Überläufen bei der Berechnung von begin und end kommen.

Ein weiteres Probleme ist, dass jeder Durchlauf den Wert von i laden muss. Dies kann je nach Architektur sehr langsam sein und somit den potentiellen Geschwindigkeitsgewinn wieder weg machen. Um dies zweifelsfrei fest zu stellen muss man die Zeit messen.

7 Akkumulation

Die Akkumulation ist die Standardanwendung der reduction-Anweisung. Es geht darum alle Elemente einer Liste auf irgendeine Weise zu kombinieren. Das Addieren aller Zahlen einer Liste ist ein Beispiel welches mit einer solchen Anweisung gelöst werden kann.

int accumulate(int arr[], int size){
	int sum = 0;
	#pragma omp parallel for reduction(+:sum)
	for(int i=0; i<size; ++i)
		sum += arr;
	return sum;
}

Leider handelt es sich jedoch nicht immer um einen Buildinoperator. Sei [i]Matrix* eine Klasse welche eine 2-dimensionale Matrix beschreibt und für welche der *= Operator entsprechend überladen wurden. Der Defaultkonstruktor initialisiert die Matrix auf die Einheitsmatrix.

class Matrix{
private:
	int 
		a,b,
		c,d
	;
public:
	Matrix():a(1), b(0), c(0), d(1){}
	Matrix&operator*=(const Matrix&right){
		Matrix left = *this;
		a = left.a*right.a + left.b*right.c;
		b = left.a*right.b + left.b*right.d;
		c = left.c*right.a + left.d*right.c;
		d = left.c*right.b + left.d*right.d;
		return *this;
	}
};

Matrix multiply(Matrix arr[], int size){
	const unsigned thread_count = omp_get_max_threads();
	Matrix*sub_prod = new Matrix[thread_count]; 
	#pragma omp parallel num_threads(thread_count)
	{
		const unsigned id = omp_get_thread_num();
		const int begin = size*id/thread_count;
		const int end = size*(id+1)/thread_count;

		for(int j=begin; j<end; ++j)
			sub_prod[id] *= arr[j];
	}
	Matrix prod;
	for(int j=0; j<thread_count; ++j)
		prod *= sub_prod[j];
	delete[]sub_prod;
	return prod;
}

Die Idee ist die Assoziativität der Multiplikation auszunutzen um das Produkt in mehrere Teilprodukte zu zerlegen und diese zuerst von jedem Thread einzeln auswerten zu lassen.

Zuerst wird bestimmt wie viele Threads überhaupt zur Verfügung stehen. Anschießend wird ein Array angelegt welches einen Eintrag für jeden Thread besitzt. In diesem Array werden die Zwischenresultate abgelegt. Die Matrizenmultiplikation ist nicht kommutativ, daher verwende ich zur Berechnung der Teilprodukte keine OpenMP-for-Schleife. Diese garantiert nämliche keine bestimmte Reihenfolge. Eine ordered OpenMP-for-Schleife ist zu steif da nicht alle Multiplikation sequentiell ablaufen müssen. Nachdem alle Teilprodukte berechnet wurden werden sie multipliziert und das Ergebnis wird zurückgegeben.

Ich habe keinen STL-vector verwendet da dieser in der Regel nicht dafür ausgelegt ist, mit mehreren Threads klar zu kommen. In diesem speziellen Fall ist es jedoch höchst unwahrscheinlich, dass dies ein Problem darstellen würde. Man könnte natürlich einen Smartpointer verwenden um die Speicherfreigabe zu garantieren. Allerdings ist das im Fall einer Ausnahme das kleinere Problem.

Die einzige***n*** Operationen die eventuell eine Ausnahme werfen könnte***n***, sind die Speicherallokierung selbst und der *= Operator. Wenn die Speicherallokierung fehlschlägt, geht kein Speicher verloren. Der *= Operator unserer Matrixklasse wirft natürlich keine Ausnahmen allerdings ist dies bei einer allgemeinen Matrixklasse nicht mehr selbstverständlich. In diesem Fall muss man allerdings sicher stellen, dass die Ausnahme sicher aus dem parallel-Block hinaus propagiert wird. Dies ist ein wesentlich größeres Problem als das kleine Speicherleck.

7.1 STL

Im Forum kommt es immer wieder zu Diskussionen zum Thema Parallelisierung und wie man die einsetzt ohne Programme von Grund auf neu zu bauen. Oft kommt dann die Idee auf die Algorithmen der STL zu parallelisieren und den Benutzer-Code unverändert zu lassen. Dies ist allerdings schwer um zu setzen. Im folgenden werde ich zeigen wie man std::accumulate parallelisieren kann und welche Probleme es dabei gibt.

Es sei folgende recht einfache Funktion gegeben:

template<class Iter, class T>
T accumulate(Iter begin, Iter end, T sum){
	for(; begin!=end; ++begin)
		sum += *begin;
	return sum;
}

Unser Ziel ist es diese Funktion durch Benutzung von OpenMP zu beschleunigen. Auf den ersten Blick erscheint dies wie eine sehr einfache Aufgabe, dies ist sie aber nicht!

Zuerst müssen wir sicherstellen, dass die zu Grunde liegenden +-Operation parallel ausgeführt werden kann. Speichert der Operator zum Beispiel Zwischenwerte in einem statischen Array, dann können wir eine Parallelisierung ohne die Operation selbst zu verändern bereits vergessen. Es ist auch möglich, dass Werte vorberechnet werden müssen und der Operator auf diese zugreift, das heißt es muss etwas initialisiert werden. Zum Beispiel eine Primzahlliste oder das Pascalsche Dreieck zur Bestimmung großer Binomialkoeffizienten sind gute Kandidaten. Dies muss vor dem ersten Aufruf passiert sein. Es darf nicht erst bis zur ersten Verwendung gewartet werden ehe man Initialisierungsarbeit erledigt.

Die +-Operation darf keine Ausnahme werfen, da C++ leider keine Möglichkeit bietet diese über Threadgrenzen hinaus zu propagieren.

Des Weiteren benötigt man die Assoziativität der +-Operation. Die Assoziativität ist notwendig damit die Summe in Teilsummen aufgespalten werden kann. Die Kommutativität ist notwendig, da nicht klar ist in welcher Reihe die Threads fertig werden. Dies erscheint offensichtlich, allerdings ist das eine zusätzliche Einschränkung gegenüber der sequentiellen Version. Oft ist sie auch nicht gegeben. Zum Beispiel ist die +-Operation von Gleitkommazahlen nicht assoziativ da die Genauigkeit von der Auswertungsreihenfolge abhängt. Normalerweise sagt man bei Gleitkommazahlen jedoch, dass die Genauigkeit nicht wichtig ist und nimmt einfach an, dass die Assoziativität gegeben ist.

Ohne die Kommutativität kann man auskommen, wenn es sein muss wie das vorherige Beispiel zeigt. Wenn sie allerdings zur Verfügung steht wird die Implementierung einfacher. Ich gehe im Folgenden davon aus, dass sie gegeben ist.

C++-Templates bieten leider praktisch keine Möglichkeiten diese Bedingungen sicher zu stellen. Die Dokumentation der Funktion ist also die einzige Stelle wo man sie angeben kann. Es bleibt nur zu hoffen, dass der Benutzer der Funktion die Dokumentation auch gründlich ließt und versteht...

Ist sicher gestellt, dass die +-Operation parallel ausgeführt werden kann müssen wir zusehen, dass wir wahlfreien Zugriff auf die Liste von Elementen haben. Dies ist nötig da es einem OpenMP-for frei steht in welcher Reihenfolge er die Schleifenduchläufe ausführen will. Ist dies nicht gegeben, kann man entweder die sequentielle Version benutzen oder die Elemente zu erst kopieren. Diese Kopie muss sequentiell erfolgen, da es keinen wahlfreien Zugriff gibt. Was besser ist hängt wieder von der +-Operation ab.

Bei Ganzzahlen ist die sequentielle Version deutlich überlegen, da die Addition wahrscheinlich so schnell ist als das Kopieren einer Zahl. Das einfache Ausrechnen ist also in etwa genauso schnell wie das Vorbereiten der parallelen Berechnung. Bei der Matrizenmultiplikation ist es allerdings besser zu kopieren. Werden Matrizen unterschiedlicher Dimensionen multipliziert so lohnt es sich unter Umständen viel mehr, zuerst dynamisch die optimale Auswertungsreihenfolge zu bestimmen als einfach nur zu parallelisieren. Ich werde im Folgenden davon ausgehen, dass es sich lohnt zu kopieren.

Der Code sieht bis hierhin in etwa folgendermaßen aus:

template<class Iter, class T>
T accumulate(Iter begin, Iter end, T init){
	return accumulate_check_category(begin, end, init, std::iterator_traits<Iter>::iterator_category());
}

template<class Iter, class T>
T accumulate_check_category(Iter begin, Iter end, T init, std::random_access_iterator_tag){
	return accumulate_parallel(begin, end-begin);
}

template<class Iter, class T>
T accumulate_check_category(Iter begin, Iter end, T init, std::random_access_iterator_tag){
	typename std::iterator_traits<Iter>::value_type T;
	std::vector<T>buffer(begin, end);
	return accumulate_parallel(buffer.begin(), buffer.size());
}

template<class Iter, class T>
T accumulate_parallel(Iter begin, int len, T init){
	// ... 
}

Die eigentliche Parallelisierung ist ein typischer Reduktionsfall. Die entsprechende OpenMP-Anweisung kann aber nur im Fall eines buildin Typen verwendet werden. Da durch die Parallelisierung Kosten entstehen lohnt es sich nicht unter einer Gewissen Grenze zu parallelisieren. Wie viele Element es geben muss, damit sich das Ganze lohnt hängt von der +-Operation, vom Prozessor und dem verwendeten Compiler ab.

Die OpenMP-Anweisung dürft bei einigen Compilern schneller sein als ihr Nachbau. Beim Nachbau muss man zwischen Nachbau mit Kommutativität und ohne Kommutativität unterscheiden. Bei der OpenMP-Anweisung braucht man dies nicht, da der buildin +-Operator kommutativ ist. Dies macht also mindestens 3 Konstanten, welche man praktisch nur empirisch bestimmen kann, die man in den Code einbauen muss. Da ich die Kommutativität angenommen habe sind es bei mir nur 2. Ich nehme einfach mal 50000 für beide an.

Ein weiteres Problem das sich stellt ist die Frage ob die Laufzeit des Operator+ von seinen Operanten abhängt oder nicht. Ob die for-Anweisung ein schedule(dynamic) benötigt oder nicht hängt davon ab. Will man alle Fälle abdecken so muss man jetzt schon 6 Konstanten bestimmen. Ich gehe mal davon aus, dass es in etwa immer gleich schnell ist. Ein Gegenbeispiel wäre die Multiplikation von Matrizen unterschiedlicher Dimensionen.

template<class Iter, class T>
T accumulate_parallel(Iter begin, int len, T init){
	accumulate_parallel_reduction(begin, len, init, boost::is_arithmetic<typename std::iterator_traits<T>::value_type>::type());
}

template<class Iter, class T>
T accumulate_parallel_reduction(Iter begin, int len, T sum, boost::true_type){
	if(len < 50000){
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}else{
		#pragma omp parallel reduction(+:ret) for
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}
	return sum;
}

template<class Iter, class T>
T accumulate_parallel_reduction(Iter begin, int len, T sum, boost::false_type){
	if(len < 50000){
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}else{
		#pragma omp parallel 
		{
			T sub_sum = T();
			for(int i=0; i<len; ++i)
				sub_sum += *(begin+i);
			#pragma omp critical
			{
				#pragma omp flush(sum)
				sum += sub_sum;
				#pragma omp flush(sum)
			}
		}
	}
	return sum;
}

Ich habe auf meinem 4-Kern-Prozessor und dem GCC-4.2 ein paar Tests gemacht um zu bestimmen wie groß bei mir diese Konstanten liegen. Unter etwa 50000 int-Elementen ist die sequentielle Version bei mir im Vorteil. Bei mehr Elementen ist die parallelisierte Version besser. Erstaunlich ist, dass die Reduktions-Answeisung in etwa genau so schnell ist wie ihr Nachbau. Ich gehe mal davon aus, dass die GCC-Entwickler den OpenMP Code noch nicht wirklich optimiert haben und in dieser Version vor allem auf korrektes Verhalten Wert gelegt haben.

8 Gibt es noch mehr?

In diesem Artikel habe ich nicht alle OpenMP-Anweisungen vorgestellt. Ich habe mich legendlich auf die mir am wichtigsten erscheinenden beschränkt. Wer OpenMP wirklich produktiv einsetzen will, dem empfehle ich die Spezifikation durch zu lesen. Für eine Spezifikation ist diese erstaunlich leicht verständlich und es sind auch viele Beispielanwendungen darin enthalten.

9 Quellen

• http://bisqwit.iki.fi/story/howto/openmp/
• http://www.openmp.org/
• http://algo2.iti.uni-karlsruhe.de/singler/mcstl/

[Ben04]

Ich bin nun wirklich kein Rechtschreib-Genie daher hab' ich auch ein paar Stellen nicht weiter beachtet (v.a. Getrennt/Zusammenschreiben von Wörtern). Ansonsten ein toller Artikel

Danke fürs Verbessern.

btw: Geht es nur mir so, dass die Smileys am Ende des Postings eingefügt werden und nicht an der Curserposition?

Ich weiß nicht was du meinst.

[Ben04]

Ich hab den Satz im flush Abschnitt umformuliert, so besser?

Es ist mir auch noch ein Fehler im Programmcode aufgefallen. Damit das Überladen auf Basis des Iteratortags funktionieren kann muss ich natürlich zwei verschiedene Tags angeben.

Inhaltsverzeichnis

• 1 OpenMP - Was ist das?
• 2 Compiler Unterstützung
• 3 Syntax
• 3.1 Parallel
• 3.2 Thread ID, single und master
• 3.3 Sections
• 3.4 Barrier
• 3.5 For
• 3.4.1 Verschachtelte Schleifen mit und ohne Collapse
• 3.4.2 Partielle sequentielle Abarbeitung mit ordered
• 3.6 Kurzschreibweise
• 4 Gemeinsamer Speicher
• 4.1 Ciritcal
• 4.2 Flush
• 4.3 Reduction
• 4.4 Locks/Mutex
• 5 Zeitmessung
• 6 Suchen in einem ungeordneten Array
• 7 Akkumulation
• 7.1 STL
• 8 Gibt es noch mehr?
• 9 Quellen

1 OpenMP - Was ist das?

OpenMP ist eine standardisierte Schnittstelle für Fortran, C und C++, welche die Entwickelung von Programmen, welche mehrere Prozessor-Kerne nutzen, vereinfachen soll. Ich werde in diesem Artikel aber nur auf die C++ Schnittstelle eingehen.

2 Compiler Unterstützung

Es gibt mehrere Versionen von OpenMP. Dieser Artikel beschäftigt sich hauptsächlich mit OpenMP 2.5. Ich werde allerdings auch auf ein paar Neuerungen in der Version 3.0 eingehen. OpenMP wird von einer Reihe von Compilern unterstützt, jedoch meistens nicht in den Basisausstattungen dieser. Mit dem GCC gibt es aber mindestens ein Compiler der jedem Leser zugänglich ist.

Folgende Tabelle liefert einen kleinen Überblick.

• GCC >= 4.2 (mit dem -fopenmp Parameter)
• MinGW Technology-Preview-Version (mit dem -fopenmp Parameter)
• VC 2005 Professional und Team System (mit dem /openmp Parameter, Standard und Express Versionen gehen nicht)

Auf der OpenMP-Webseite findet man eine detailliertere Tabelle.

3 Syntax

OpenMP besteht zu einem Großteil aus Compiler-Erweiterungen. Es gibt eine Reihe von Pragma-Anweisungen welche die OpenMP-Funktionalität bereitstellen. Jede von ihnen hat folgenden Syntax:

#pragma omp Anweisung
C++ Code auf den sich die Anweisung bezieht

Neben diesen pragma-Konstrukten gibt es noch eine Reihe von normalen Funktionen. Diese werden im Header omp.h definiert. Alle haben den Präfix omp_.

Jeder Compiler der sich an den C++-Standard hält muss unbekannte Pragma-Anweisungen ignorieren und höchstens mit einer Warnung versehenen. Fehler sich nicht standard-konform. Durch das Benutzen eines Dummy-Headers, ist es deswegen in vielen Fällen möglich eine OpenMP-Anwendung sequentiell auf nicht OpenMP-fähigen Compilern zu übersetzen. Dies ist auch empfehlenswert wenn der Compiler zwar OpenMP unterstützt, allerdings der Zielprozessor keine echte Parallelisierung unterstützt.

In diesem Artikel werde ich versuchen die Anweisungen welche kein sequentielles Fallback-Konstrukt haben zu ignorieren. Sie bieten auch eher selten Funktionalität an, welche man nicht auch anders erhalten kann. Daher ist dies an sich kein großer Verlust.

3.1 Parallel

Die grundlegend Anweisung, ohne die bei OpenMP nichts auskommt, ist die parallel-Anweisung. Sie erzeugt eine Reihe von Threads und führt den nachfolgenden C++ Block mit jedem Thread aus. Dies bedeutet dass der gleiche Code mehrmals gleichzeitig ausgeführt werden kann und zwar von jedem Thread einmal. Wie viele Threads erstellt werden, ist der Implementierung von OpenMP freigestellt. Meistens werden aber so viele Threads erzeugt wie es im Rechner Kerne gibt.

Es ist also in der Regel nicht möglich vorher zu sagen wie oft ein parallel-Block ausgeführt wird.

int main(){
	printf("Hier ist noch alles einfach und sequentiell.\n");
	#pragma omp parallel
	{
		printf("Hallo ");
		printf("Welt!\n");
	}
	printf("Alle Threads haben Hallo gesagt.");
}

Im parallel-Block ist ganz normaler C++-Code erlaubt. Man kann C++-Objekte anlegen und auch Funktionen aufrufen. Sogar Ausnahmen sind erlaubt, solange sie nicht aus dem parallel-Block rausfliegen.

Solange sich die Threads nicht in die Quere kommen, ist das auch schon alles was man machen muss um sämtliche Kerne seines Rechners zu benutzen. Allerdings ist dies leider nur sehr selten der Fall. Bereits in diesem einfachen Beispiel ist es möglich, dass Threads sich streiten. Wenn zwei Threads gleichzeitig etwas ausgeben, ist es nämlich nicht klar welcher zuerst dran kommt. Alles von Funktionieren-wie-erwartet bis hin zu Abstürzen ist eine möglich Konsequenz.

Was das Problem noch komplizierter macht, ist dass das Verhalten nicht vollständig deterministisch ist. Wie viele Threads erzeugt werden kann von Lauf zu Lauf unterschiedlich sein und welcher Thread wann dran kommt kann man noch schlechter vorhersagen. Dies führt dazu, dass manche Bugs nur mit einer gewissen Wahrscheinlichkeit eintreten. Das Debuggen solcher Fehler ist sehr schwer.

An sich müsste ich hier sicherstellen, dass zu jedem Zeitpunkt immer höchstens nur ein Thread etwas auf die Konsole schreibt. Bei printf kann man aber normalerweise davon ausgehen, dass legendlich die Reihenfolge der Ausgaben etwas durcheinander geraten kann. Abstürze gehören eher in den Bereich der Theorie. printf leidet normalerweise weniger unter diesem Problem als std::cout, deswegen werde ich in diesem Artikel ausnahmsweise printf bevorzugen. Im Folgenden werde ich davon ausgehen, dass Ausgaben mit printf einfach funktionieren. Alles andere würde die Beispiele zu kompliziert machen.

Da man nicht wirklich weiß wie oft ein parallel-Block ausgeführt wird, ist eine Parallelisierung recht umständlich ohne weitere OpenMP-Konstrukte. Er eignet sich aber in Regel vorzüglich um threadspezifische Variablen zu initialisieren.

Mit Hilfe der num_threads(n)-Anweisung kann man festlegen wie viele Threads erstellt werden. n muss eine Ganzzahl sein und braucht nicht zum Zeitpunkt der Übersetzung bekannt zu sein.

int main(){
	printf("Wie viele Threads sollen es sein?\n");
	int n;
	scanf("%d", &n);
	#pragma omp parallel num_threads(n)
	{
		printf("Hallo ");
		printf("Welt!\n");
	}
	printf("Alle Threads haben Hallo gesagt.");
}

Wenn man für n eine große Zahl eingibt, kann man mit großer Wahrscheinlichkeit beobachten wie sich die Threads in die Quere kommen. Es folgt eine Beispielausgabe, so wie ich sie beobachten konnte.

Hallo Hallo Welt!
Hallo Welt!
Welt!
Hallo Welt!

Während der eine Thread noch am Schreiben ist, funkt ein anderer ihm einfach dazwischen. Versuchen Sie nach zu vollziehen wie diese Ausgabe zustande kam. Es ist sehr wichtig dieses Problem zu verstehen, denn sonst werden Sie mit OpenMP oder auch Multithread-Programmierung im Allgemeinen nur Ärger haben.

omp_get_max_threads gibt an wie viele Threads im nächsten parallel-Block erzeugt werden.

Bei den meisten OpenMP-Implementierungen werden die Threads nur einmal erstellt. Dies gilt auch wenn es mehrere parallel-Blocks gibt. Zwischen den Blocks gehen die zusätzlichen Threads in einen Wartezustand. Sie werden erst am Ende des Programms zerstört. Die in Regel recht hohen Threaderstellungskosten muss man also nur genau einmal zahlen.

parallel-Blöcke dürfen verschachtelt werden. In der Standardeinstellung erzeugt der verschachtelte Block aber keine neuen Threads.

3.2 Thread ID, single und master

Jeder Thread besitzt eine eindeutige Identifikationsnummer. Es handelt sich dabei um fortlaufende Nummern und nicht etwa um gecastete Zeiger oder ähnliches. Man kann die ID also als Index in einem Array benutzen. Den Thread mit der ID 0 nennt man auch noch Master.

Mit Hilfe von omp_get_thread_num kann man die ID eines Threads bestimmen und durch einen Aufruf von omp_get_num_threads kann man die Gesamtanzahl der Threads, welche vom parallel-Block erzeugt wurden, heraus finden.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());	
		if(omp_get_thread_num() == 0)
			printf("Es gibt %d Threads.\n", omp_get_num_threads());		
	}

}

Wie man an diesem Beispiel sieht, ist es möglich auch in einem parallel-Block, Instruktionen nur einmal aus zu führen. Allerdings gibt es dafür bessere OpenMP-Konstrukte.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());	

		#pragma omp master
		{
			printf("Es gibt %d Threads.\n", omp_get_num_threads());	
		}

		#pragma omp single
		{
			printf("Thread %d war hier.\n", omp_get_thread_num());	
		}

		#pragma omp single
		{
			printf("Danach war Thread %d da.\n", omp_get_thread_num());	
		}
	}
}

Ein master-Block wird nur vom Master-Thread ausgeführt. Der Code ist äquivalent zur if-Abfrage aus dem vorherigen Beispiel. Ein single-Block wird genau einmal ausgeführt, allerdings ist es egal von welchem Thread.

Sämtliche Threads warten am Ende eines single-Blocks bis auch der letzte es bis dahin geschafft hat. Dadurch ist sichergestellt, dass der zweite single Block nicht ausgeführt wird ehe der erste vollständig abgearbeitet wurde. Durch eine nowait Anweisung kann man das Warten unterdrücken.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			printf("Manchmal bin ich als erstes da.\n");	
		}

		#pragma omp single nowait
		{
			printf("Hin und wieder gewinne aber auch ich das Rennen.\n");	
		}
	}
}

Wenn man nowait verwendet dann wird nur noch garantiert, dass die Ausführung des ersten Blocks vor der des zweiten beginnt. Es ist aber gut möglich, dass der zweite vor dem ersten beendet wird. Es ist auch möglich, dass die Ausführung des ersten Threads noch vor der ersten Anweisung im Block unterbrochen wird und es dadurch erscheint, als wenn der zweite Block vor dem ersten gestartet werden würde. In der Praxis kann man also nicht darauf bauen, dass der erste vor dem zweiten gestartet wird, auch wenn dies mit einer hohen Wahrscheinlichkeit so sein wird.

Am Ende eines Master-Block wird nicht gewartet.

Bei verschachtelten parallel-Blöcken werden die IDs durch den innersten Block neu verteilt. In der Standardeinstellung erhalten alle Threads die ID 0 und die Anzahl der Threads wird als 1 gemeldet. Die Idee dahinter ist, dass sich alles auf den innersten Block bezieht. Dieser hat aber keinen Thread erzeugt und folglich gibt es nur den der den Block angestoßen hat. Eine weitere Konsequenz ist, dass single- und master-Blöcke mehrfach ausgeführt werden. Ein verschachtelter parallel-Block wird demnach also in der Standardeinstellung sequentiell abgearbeitet. Mit num_threads(n) kann das erzeugen weiterer Threads erzwingen.

Die Problematik wird durch folgendes Beispiel gut illustriert.

int main(){
	#pragma omp parallel
	{
		#pragma omp parallel
		{
			printf("omp_get_thread_num() = %d,"
				" omp_get_num_threads() = %d\n", 
				omp_get_thread_num(), 
				omp_get_num_threads());	
		}
	}
}

3.3 Sections

Alternativ zu den single-Blöcken kann man auch sections verwenden.

int main(){
	#pragma omp parallel
	{
		#pragma omp sections
		{
			{
				printf("Thread %d war hier.\n", omp_get_thread_num());			
			}
			#pragma omp section
			{
				printf("Aber hier war Thread %d.\n", omp_get_thread_num());		
			}
			#pragma omp section
			{
				printf("und eine weitere section");		
			}
		}
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

Jede der drei Sektionen wird genau einmal ausgeführt. Am Ende des sections-Block warten die Threads ehe sie fortfahren. Auch hier kann man das Warten durch nowait unterdrücken.

Man könnte obiges Beispiel auch mit single ausdrücken.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			printf("Thread %d war hier.\n", omp_get_thread_num());			
		}
		#pragma omp single nowait
		{
			printf("Aber hier war Thread %d.\n", omp_get_thread_num());		
		}
		#pragma omp single nowait
		{
			printf("und eine weitere section");		
		}
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

Bei der single-Variante ist der erste Block mit hoher Wahrscheinlichkeit der erste der gestartet wird. Bei der sections-Variante sind alle Blöcke gleich wahrscheinlich.

3.4 Barrier

Wenn an einer beliebigen Stelle gewartet werden soll so kann man eine barrier-Anweisung benutzen.

int main(){
	#pragma omp parallel
	{
		printf("Thread %d sagt Hallo \n", omp_get_thread_num());
		#pragma omp barrier
		printf("Thread %d sagt Tschüss \n", omp_get_thread_num());
	}
}

In diesem Beispiel ist sicher gestellt, dass alle Threads "Hallo" sagen ehe sich der erste verabschiedet.

3.5 For

Mit Hilfe von OpenMP kann man auch einfache Schleifen parallelisieren.

int main(){
	#pragma omp parallel
	{
		#pragma omp for
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);
	}
}

Anders als bei einer normalen Schleife, werden die Durchläufe hier nicht unbedingt nacheinander abgearbeitet, sondern möglicherweise gleichzeitig und auch in zufälliger Reihenfolge. Dies führt dazu, dass die Ausgabe des obigen Codes nicht immer die gleiche ist. Die Reihenfolge der Zeilen kann beliebig permutiert sein.

Um dies zu bewerkstelligen werden die Durchläufe möglichst gleichmäßig zwischen den Threads verteilt. Zum Beispiel bei zwei Threads, könnte der erste die Durchläufe mit 0<=i<=4 abarbeiten und der andere die restlichen mit 5<=i<=9. Die genaue Verteilung ist dem Compiler frei gestellt. Sie hängt jedoch nur von der Anzahl der Threads und der Anzahl der Schleifendurchläufen ab, das heißt wenn aus irgendeinem Grund ein Thread massive Verspätung hat, dann wird nicht kurzerhand die Arbeit umgeschichtet.

int main(){
	#pragma omp parallel
	{
		// Sorgt für massive Verspätung eines Threads
		#pragma omp single nowait
		{
			for(int i=0; i<INT_MAX; ++i)
				for(int j=0; j<INT_MAX; ++j)
					{}
		}

		// Tja hier warten dann die anderen...
		#pragma omp for
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);
	}
}

Obwohl der eine Thread sehr lange mit dem single-Block beschäftigt ist, kommt keiner der wartenden Threads auf die Idee bei der for-Schleife Überstunden zu leisten. Dies ist nur mit einer dynamischen Arbeitsverteilung möglich.

int main(){
	#pragma omp parallel
	{
		#pragma omp single nowait
		{
			for(int i=0; i<INT_MAX; ++i)
				for(int j=0; j<INT_MAX; ++j)
					{}
		}

		#pragma omp for schedule(dynamic)
		for(int i=0; i<10; ++i)
			printf("%d^2 = %d\n", i, i*i);

		// Hier wird erst gewartet wenn die Schleife vollständig abgearbeitet wurde.
	}
}

Wenn ein Thread fertig ist, schaut er ob es noch Arbeit gibt und übernimmt diese dann auch gleich. Ein ähnliches Problem ergibt sich wenn die Durchläufe unterschiedlich lang dauern.

Am Ende eines for-Blocks wird gewartet. Dies kann man wie üblich mit nowait unterdrücken.

Da schedule(dynamic) einen größere Verwaltungskosten als die normale Version hat, kann es sein, dass die normale Version in manchen Fällen schneller ist, obwohl nicht alle Threads immer beschäftigt sind. Die normale Version kommt nämlich ganz ohne Kontakt zwischen den einzelnen Threads zurecht und ist daher sehr schnell. Ich werde weiter unten auf konkrete Implementierungsmöglichkeiten eingehen.

Damit eine Schleife überhaupt parallelisiert werden kann muss sie eine Reihe von sehr strikten Kriterien erfüllen.

• Die Schleifenvariable muss eine Ganzzahl sein. In OpenMP 3.0 sind auch Random-Access-Iteratoren und Zeiger erlaubt.
• Die Abbruchbedingung muss ein Ordnungsvergleich (<,>, <= oder sein und auf der linken Seite muss sich die Schleifenvariable und nur die befinden. Die rechte Seite darf nicht während der Ausführung der Schleife verändert werden.
• Die Inkrementierung muss in immer gleich großen Schritten stattfinden. Dies bedeutet, dass man nur --, ++, += oder -= verwenden darf um die Schleifenvariable zu verändern. Die Größe der Schritte darf nicht während der Ausführung der Schleife verändert werden.
• break und return sind verboten. continue ist erlaubt. goto darf nur verwendet werden wenn man nicht aus der Schleife raus springt.

Folgende Schleifen lassen sich nicht mit OpenMP parallelisieren.

int n = foo();

// Auf der einen Seite der Bedingung muss ein einfaches i stehen.
for(int i=0; i*i<n; ++i)
	printf("%d", i);

// Die Inkrementierungsschritte sind nicht gleichgroß
for(int i=0; i<n; i*=3)
	printf("%d", i);

// Dies ist erlaubt, allerdings ist unspezifiziert wie oft bar und foobar aufgerufen werden.
for(int i=-10; i<n*n*bar(); i+=foobar(n))
	printf("%d", i);

// Man darf nicht aus der Schleife raus springen.
int a[10] = { ... }
for(int j=0; j<10; ++j)
	if(a[j] == n)
		return i;

// n darf nicht in der Scheife verändert werden.
for(int i=0; i<n*10; i+=n)
	++n;

3.4.1 Verschachtelte Schleifen mit und ohne Collapse

OpenMP-for-Schleifen dürfen nicht verschachtelt werden. Man kann allerdings mit einem zusätzlichen parallel-Block Abhilfe schaffen. Dies ist jedoch selten sinnvoll da die Schleife zwar eine OpenMP-Schleife ist, aber dennoch sequentiell abgearbeitet wird.

#pragma omp parallel
{
	#pragma omp for
	for(int i=0; i<10; ++i)
	{
		#pragma omp parallel
		{
			#pragma omp for
			for(int j=0; j<10; ++j)
			{
				printf("%d*%d = %d\n", i, j, i*j);
			}
		}
	}
}

Es wird aber nur die äußere Schleife parallelisiert und die innere sorgt nur für Probleme und Kosten. Besser ist:

#pragma omp parallel
{
	#pragma omp for
	for(int i=0; i<10; ++i)
	{
		// Dieses for ist ein ganz normales C/C++ for
		for(int j=0; j<10; ++j)
		{
			printf("%d*%d = %d\n", i, j, i*j);
		}

	}
}

In OpenMP 3.0 gibt es eine spezielle collapse-Schleife um über ein mehrdimensionales Array zu iterieren. Dies entspricht einem der häufigsten Einatzgebiet von mehrfach verschachtelten Schleifen. Die Schleifen müssen allerdings sehr einfach gehalten sein, damit eine Parallelisierung überhaupt erfolgen kann. Alle Schleifenvariablen müssen völlig unabhängig von einander sein. Es darf nur die innerste Schleife beliebigen Code enthalten. Alle anderen dürfen nur verschachtelte Schleifen enthalten. Das Argument der collapse-Anweisung gibt an wie tief die Schleifen verschachtelt sind.

Dies sind sehr starke Einschränkungen und deswegen deckt das folgende Beispiel an sich auch quasi jede erlaubte Variante einer collapse-Schleife ab.

int main(){
	#pragma omp parallel
	{
		int foo[10][3][7];
		#pragma omp for collapse(3)
		for(int l=0; l<10; ++l)
			for(int j=0; j<3; ++j)
				for(int k=0; k<7; ++k){
					// Nur hier darf beliebiger Code stehen!
					foo[l][j][k] = 0;
				}

}

3.4.2 Partielle sequentielle Abarbeitung mit ordered

Man kann in eine OpenMP-for Schleife einen ordered-Block packen. Dieser sorgt dafür, dass ein Teil der Schleife sequentiell abgearbeitet wird.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<=10; ++i){
			double v = i;
			for(int j=1; j<10; ++j)
				v *= i;
			#pragma omp ordered
			{
				printf("%d^10 = %.0f\n", i, v);
			}
		}
	}
}

Die Berechnung der Werte kann in beliebiger Reihenfolge und parallel erfolgen. Die Ausgabe ist jedoch sequentiell, das heißt die Ausgabe ist immer die gleiche.

0^10 = 0
1^10 = 1
2^10 = 1024
3^10 = 59049
4^10 = 1048576
5^10 = 9765625
6^10 = 60466176
7^10 = 282475249
8^10 = 1073741824
9^10 = 3486784401
10^10 = 10000000000

Ein ordered-Block darf nur in einer speziellen ordered-for-Schleife auftreten. Jeder Durchlauf darf nur maximal einem ordered-Block begegnen. Folgendes Beispiel ist demnach falsch.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<10; ++i){
			#pragma omp ordered
			printf("%d^10 =", i);

			double v = i;
			for(int j=1; j<10; ++j)
				v *= i;

			// Hier wird etwas schief gehen da 
			// dies der zweite Block ist.
			#pragma omp ordered
			printf("%.0f\n", v);
		}
	}
}

Dies bedeutet jedoch nicht, dass man nur einen ordered-Block pro Schleife haben darf.

int main(){
	#pragma omp parallel
	{
		#pragma omp for ordered
		for(int i=0; i<10; ++i){
			if(i%3 == 0){
				#pragma omp ordered
				printf("%d ist ein Vielfaches von 3\n", i);
			}else{
				#pragma omp ordered
				printf("%d ist kein ein Vielfaches von 3\n", i);
			}
		}
	}
}

In dem Beispiel ist sichergestellt, dass kein Durchlauf zwei verschiedene Blöcke zu sehen bekommt.

3.6 Kurzschreibweise

Wenn ein parallel-Block nur eine for-Schleife enthält, dann gibt es eine Kurzschreibweise.

#pragma omp parallel for
for(int i=0; i<10; ++i)
{
	for(int j=0; j<10; ++j)
	{
		printf("%d*%d = %d\n", i, j, i*j);
	}

}

Analog gibt es auch ein parallel sections-Block.

Folgendes Beispiel berechnet alle Primzahlen kleiner als eine gegebene Zahl.

bool is_prime(unsigned n){
	for(unsigned i=2; i<n; ++i)
		if(n%i == 0)
			return false;
	return true;
}

int main(){
	unsigned n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%u", &n);

	// schedule(dynamic) ist notwendig da is_prime keine
	// konstante Laufzeit besitzt. 
	#pragma omp parallel for schedule(dynamic)
	for(int i=2; i<=n; ++i)
		if(is_prime(i))
			printf("%u ist prim\n", i);
}

4 Gemeinsamer Speicher

Die Zugriffsregeln für Variablen sind die die man erwarten würde. Variablen leben so lange wie der Block in dem sie definiert wurden. Sie können verdeckt werden und man kann auf Variablen einer höheren Ebene zugreifen.

int a;

int main(){
	int b;
	#pragma omp parallel
	{
		int c;
		// a und b werden geteilt, jeder Thread hat sein eigenens c
	}
}

Wenn ein Thread schreibend auf a oder b zugreift dann ergeben sich sämtliche vom Multithreading her bekannte Probleme. OpenMP bietet aber leider nur wenig Revolutionäres um dieses alt bekannte Problem zu lösen.

4.1 Ciritcal

Ein critical-Block ist ein Block in dem sich immer nur genau ein Thread befindet. Wenn bereits ein Thread im Block ist, wartet jeder Thread der hinein will, bis der Block wieder frei wird.

Betrachten Sie folgendes Beispiel:

int main(){
	#pragma omp parallel for
	for(int i=0; i<10; ++i){
		int n;
		#pragma omp critical
		{
			printf("Geben sie eine Zahl ein: ");
			scanf("%d", &n);
		}
		printf("%d * %d = %d\n", i, n, i*n);
	}
}

Hier wird ein critical-Block verwendet um sicher zu stellen, dass der Benutzer nie aufgefordert wird zwei oder mehr Zahlen gleichzeitig einzugeben. Nimmt man das critical weg dann kann es zu folgender Ausgabe kommen.

Geben sie eine Zahl ein: Geben sie eine Zahl ein:

Jeder critical-Block besitzt einen Namen. Wenn zwei Blöcke denselben Namen besitzen, kann nur genau ein Thread in einem der beiden sein. Zwei Blöcke mit verschiedenen Namen beeinflussen sich nicht. Das heißt es ist durchaus möglich, dass 2 Threads in critical-Blöcken sind sofern sie verschiedene Namen haben.

Wenn bei einem critical-Block kein Name angegeben wird, so wie im ersten Beispiel, besteht der Name aus dem leeren Wort. Alle leeren Wörter werden als gleich angesehen. Im Klartext bedeutet dies, dass keine zwei Threads gleichzeitig in irgendeinem namenlosen Block sein können.

Leider sind die Namen global und liegen in keinem Namensraum. Dadurch sind Namenskonflikte möglich. Es ist allerdings sehr unwahrscheinlich, dass dies andere Auswirkungen hat als einfach nur ein langsameres Programm, da lediglich ein Thread warten muss obwohl dies nicht nötig gewesen wäre. Es ist sichergestellt, dass sich normalen Symbolen wie Funktionsnamen nicht in die Quere kommen.

int main(){
	FILE*file = fopen("out.txt", "w");
	#pragma omp parallel for
	for(int i=0; i<10; ++i){
		int n;

		#pragma omp critical(console)
		{
			printf("Geben sie eine Zahl ein : ");
			scanf("%d", &n);
		}

		#pragma omp critical(console)
		{
			printf("%d * %d = %d\n", n, n, n*n);
		}

		#pragma omp critical(file)
		{
			fprintf(file, "%d * %d = %d\n", n, n, n*n);
			fprintf(file, "%d + %d = %d\n", n, n, n+n);
		}

		#pragma omp critical(console)
		{
			printf("%d + %d = %d\n", n, n, n+n);
		}
	}
	fclose(file);
}

In diesem Beispiel ist sicher gestellt, dass nichts auf die Konsole geschrieben wird, falls der Benutzer aufgefordert wurde etwas einzugeben. In der Datei müssen beide zueinander passenden Rechnungen hinter einander ausgegeben werden. Auf der Konsole ist dies nicht garantiert.

4.2 Flush

Viele Rechner besitzen Register und mehrere Level von Speicher-Caches, um den Zugriff auf Variablen zu beschleunigen. Anstatt direkt mit der Variable im Hauptspeicher zu arbeiten wird diese zuerst in einen schnelleren Cache-Speicher kopiert. Anschließend wird mit der Kopie gerechnet und nach einiger Zeit wird die Variable in den Hauptspeicher zurück kopiert. Dies führt dazu, dass es mehrere Versionen von der gleichen Variable gibt welche nicht immer gleich sind.

In einem sequentiellen Programm stellt dies kein Problem dar, da zu jedem Zeitpunkt bekannt ist wo sich die aktuelleste Version der Variable befindet. Jedoch ist dies bei mehreren Threads nicht mehr so einfach. Der eine weiß nicht wann der andere fertig ist und wann er die Werte synchronisieren muss.

Um dieses Problem zu lösen bietet OpenMP die flush-Anweisung an. Sie synchronisiert die lokale Version des Threads mit der des Hauptspeichers. Dies bedeutet, dass jeder Thread nachdem er fertig ist, eine geteilte Variable zu verändern, diese auch flushen muss. Ein Thread der lesend auf eine Variable zugreift will, muss die Variable auch flushen damit er sicher ist die neuste Version zu Gesicht zu bekommen. Wenn ein Thread auch mit einer veralteten Version arbeiten kann, braucht er die Variable nicht zu flushen. Wenn sicher ist, dass eine Variable in keinem parallel-Block verändert wird, dann kann man sich das flush auch sparen. Es ist nämlich egal wann ein Thread seine lokale Version aus dem Hauptspeicher lädt da er immer die gleichen Daten lesen würde.

Folgendes Programm berechnet das Faktoriell einer Zahl.

int main(){
	unsigned n, fact = 1;
	printf("Geben sie eine Zahl ein: ");
	scanf("%u", &n);

	#pragma omp parallel
	{
		unsigned local_fact = 1;

		// Auf n kann man ohne flush zugreifen,
		// da es nicht verändert wird.
		#pragma omp for
		for(int i=1; i<=n; ++i){
			// local_fact braucht nicht geflusht 
			// zu werden, da es nicht geteilt wird.
			local_fact *= i;
		}

		#pragma omp critical
		{
			// Zuerst holt der Thread sich die neuste Version
			// aus dem Hauptspeicher.
			#pragma omp flush(fact)

			// Danach wird die Variable verändert.
			fact *= local_fact;

			// Anschließend wird die Variable wieder in den
			// Hauptspeicher geschrieben, damit andere Threads
			// sich die neuste Version von dort laden können.
			#pragma omp flush(fact)
		}	
	}

	printf("%u! = %u\n", n, fact);
}

Ein fehlendes flush bedeutet nicht, dass die Variablen nicht synchronisiert werden. Es steht der Implementierung frei dies zu tun wann immer sie will. Dies führt leider oft zu nur sehr schwer reproduzierbaren Bugs. Zum Beispiel ist es sehr wahrscheinlich, dass das Betriebssystem die Rechnerzeit anders auf die Threads verteilt, wenn die Hardware unter voller Last läuft, als wenn es nur mit kleinen Texteingaben konfrontiert wird. Es ist daher sinnvoll im Zweifelsfall lieber ein flush zuviel zu schreiben.

4.3 Reduction

Mit reduction lassen sich Berechnung vom Typ s += f(1)+f(2)+...+f(n) parallelisieren. Jeder Thread rechnet ein paar Terme aus und summiert die auf welche er berechnet hat. Nachdem alle Threads fertig sind zählt jeder Thread sein Zwischenergebnis zur Variable s hinzu.

Diese Zwischenergebnis wird in einer Variable mit dem gleichen Namen wie s gespeichert. Diese Hilfsvariable wird im reduction-Block angelegt und mit dem neutralen Element (0 im Fall der Addition) initialisiert. Man sollte sich aber nicht auf dieses Verhalten verlassen. Wenn ein Programm sequentiell übersetzt wird, das heißt ohne OpenMP Unterstützung, dann gibt es keine Zwischenergebnisse und deswegen hat s beim ersten Durchlauf auch nicht unbedingt den Wert des neutralen Elements.

Neben + lassen sich auch -, *, &, |, ^, && und || auf diese Weise parallelisieren. Es muss bei den Operatoren jedoch um ihre Buildin-Version handeln.

• +, -, | und ^ besitzen das neutrale Element 0.
• * besitzt das neutrale Element 1.
• & besitzt das neutrale Element ~0. (Bei ~0 sind alle Bits gesetzt sind.)
• && besitzt das neutrale Element true.
• || besitzt das neutrale Element false.

Bei Gleitkommazahlen hängt die Genauigkeit von der Auswertungsreihenfolge ab. Die Reihenfolge ist nicht spezifisiert, also auch nicht die Genauigkeit des Ergebnisses. Es ist nicht einmal garantiert, dass es immer das gleiche Ergebnis sein wird. Gleitkommazahlen in Kombination mit einer Reduktion haben einen gewissen nicht deterministischen Teil im Ergebnis.

In der Fortran Version von OpenMP gibt es zusätzlich noch Min- und Max-Reduktionen. Die wurden aus mir unerklärlichen Gründen in der C/C++ Version vergessen.

Folgendes Beispiel bestimmt ob eine Zahl prim ist.

int main(){
	int n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%d", &n);

	// Der Initialwert von is_prime ist wichtig da dieser auch
	// in die Berechnung einbezogen wird.
	bool is_prime = true;

	#pragma omp parallel reduction(&&: is_prime)
	{
		// jeder bekommt sein eigenes is_prime,
		// welches mit true, dem neutralen Element von &&,
		// initialisiert wurde.

		#pragma omp for
		for(int i=2; i<n; ++i)
			if(n%i == 0)
				is_prime = false;

		// Alle is_prime (auch das aus der main Funktion) werden 
		// mittels && verknüpft. Im Klartext bedeutet dies, dass 
		// falls eines false ist, ist das Resultat false.
	}

	if(is_prime)
		printf("%d ist prim\n", n);
	else
		printf("%d ist nicht prim\n", n);

}

Am Ende des parallel Blocks wird folgende Rechnung ausgeführt:

is_prime = is_prime && thread[0].is_prime && thread[1].is_prime && ...  && thread[thread_num-1].is_prime

Bei der thread[j]-Notation handelt es sich legendlich um Pseudocode.

Auch hier gibt es eine Kurzschreibweise. Folgendes Beispiel berechnet das Faktoriell einer Zahl.

int main(){
	unsigned n;

	printf("Geben sie eine Ganzzahl ein: ");
	scanf("%u", &n);

	unsigned fact = 1;

	#pragma omp parallel for reduction(*: fact)
	for(int i=2; i<n; ++i)
		fact *= i;

	printf("%u! = %u\n", n, fact);
}

reduction erlaubt es dem Compiler eine Reihe von Optimierungen durchzuführen, welche bei der äquivalenten Formulierungen mit critical und flush nur sehr schwer möglich sind.

4.4 Locks/Mutex

OpenMP bietet auch Mutexe an. Hier weicht OpenMP an sich nicht von allgemein bekanten Konzepten ab. Es gibt einen Lock-Type namens omp_lock_t. Jeder Mutex ist entweder offen oder geschlossen. Folgende Funktionen werden bereitgestellt um damit zu arbeiten:

• void omp_init_lock(omp_lock_t*)
  Initialisiert den Mutex. Der Mutex ist danach offen.
• void omp_destroy_lock(omp_lock_t*)
  Zerstört den Mutex.
• void omp_set_lock(omp_lock_t*)
  Versucht den Mutex zu schließen. Wurde der Mutex bereits von einem anderen Thread geschlossen so wird gewartet bis, dass dieser ihn wieder geöffnet hat.
• void omp_unset_lock(omp_lock_t*)
  Öffnet den Mutex wieder
• int omp_test_lock(omp_lock_t*)
  Versucht den Mutex zu schließen. Wurde der Mutex bereits von einem anderen Thread geschlossen so wird 0 zurück gegeben, ansonsten 1.

omp_lock_t ist ein nicht rekursives Mutex. Das heißt ein Thread darf einen Mutex nicht mehrmals schließen. Er würde in dem Fall nämlich auf sich selbst warten. Es gibt aber auch eine rekursive Variante und die heißt omp_nest_lock_t. Die Funktionen um damit umzugehen heißen ähnlich und funktionieren analog.

• void omp_init_nest_lock(omp_next_lock_t*)
• void omp_destroy_nest_lock(omp_next_lock_t*)
• void omp_set_nest_lock(omp_next_lock_t*)
• void omp_unset_nest_lock(omp_next_lock_t*)
• int omp_test_nest_lock(omp_next_lock_t*)

Da diese Schnittstelle alles anderes als modernes C++ ist, bietet es sich an diese Funktionen zu wrappen.

class Mutex
{
public:
	Mutex()				{omp_init_lock(&lock);}
	~Mutex()			{omp_destroy_lock(&lock);}
	void lock()			{omp_set_lock(&lock);}
	void unlock()			{omp_unset_lock(&lock);}
	bool try_to_lock()		{return omp_test_lock(&lock);}
private:
	Mutex(const Mutex&);
	Mutex&operator=(const Mutex&);
	mp_lock_t lock;
};

class ScopedLock
{
public:
	ScopedLock(Mutex&lock)		:lock(lock){lock.lock();}
	~ScopedLock()			{lock.unlock();}
private:
	ScopedLock(const ScopedLock&);
	ScopedLock&operator=(const ScopedLock&);
	Mutex&lock;
};

5 Zeitmessung

Neben den Funktionen welche direkt mit der Multithread-Programmierung bietet OpenMP noch zwei nützliche Funktionen zur Zeitmessung an.

• double omp_get_wtime();
  Gibt die momentane Zeit in Sekunden zurück. Es ist nicht definiert wann der Zeitpunkt 0 ist. Es ist legendlich bekannt, dass es nicht während dem Lauf eines Programms verändert wird und in der Vergangenheit liegt. Dies ist aber nicht weiter schlimm, wenn man die Funktion nur für Zeitmessungen einsetzen will. Die Genauigkeit dieser Funktion ist nicht spezifisiert und darf leicht von Thread zu Thread abweichen. Man kann die Genauigkeit aber mit omp_get_wtick in Erfahrung bringen.
• double omp_get_wtick();
  Gibt den kleinstmöglichen Abstand zwischen 2 verschiedenen von omp_get_wtime zurückgegebenen Zeitpunkten an. Der Rückgabewert von omp_get_wtick darf von Thread zu Thread variieren.

6 Suchen in einem ungeordneten Array

Das Suchen in einem unsortierten Array lässt sich sehr leicht mit OpenMP beschleunigen.

// Ohne OpenMP
int find_first(int arr[], int size, int value){
	for(int j=0; j<size; ++j)
		if(arr[j] == value)
			return j;
	return -1;
}

// Mit OpenMP
int find_any(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel for
	for(int j=0; j<size; ++j)
		if(arr[j] == value){
			i = j;
			#pragma omp flush(i)		
		}
	return i;
}

Das aktuelle i braucht nicht zuerst per flush-Anweisung geladen zu werden, da der Wert eh nur überschrieben werden würde.

Wie man bereits am Namen erkennt, ist bei der parallelen Version nicht definiert welchen Wert gefunden wird, wenn es mehrere gleiche Werte im Array gibt. Wenn dies jedoch nötig ist, kann man dieses Problem leicht mit ordered beheben.

int find_first(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel for ordered
	for(int j=size-1; j>=0; --j)
		if(arr[j] == value){
			#pragma omp ordered
			i = j;
			#pragma omp flush(i)		
		}
	return i;
}

Eines der Nachteile der parallelen Version ist, dass immer alle Elemente angeschaut werden. Dies ist aber nicht immer notwendig denn man kann die Schleife abbrechen sobald ein Element gefunden wurde. OpenMP bietet allerdings keine Möglichkeit dies zu tun. Man kann versuchen die Schleife nach zu bauen.

int find_any(int arr[], int size, int value){
	int i = -1;
	#pragma omp parallel
	{
		const int begin = size*omp_get_thread_num() / omp_get_num_threads();
		const int end = size*(omp_get_thread_num()+1) / omp_get_num_threads();

		for(int j=begin; j<end; ++j){
			if(arr[j] == value)
				i = j;
			#pragma omp flush(i)
			if(i != -1)
				break;
		}
	}
	return i;
}

Dieser Code dürfte recht nahe an dem sein, was der Compiler aus einer OpenMP-Schleife mit statischer Verteilung macht. In der Regel sollte man aber die OpenMP-for-Schleife bevorzugen. Diese ist weit weniger fehleranfällig. So kann es beim Beispiel oben zum Beispiel zu Überläufen bei der Berechnung von begin und end kommen.

Ein weiteres Probleme ist, dass jeder Durchlauf den Wert von i laden muss. Dies kann je nach Architektur sehr langsam sein und somit den potentiellen Geschwindigkeitsgewinn wieder weg machen. Um dies zweifelsfrei fest zu stellen muss man die Zeit messen.

7 Akkumulation

Die Akkumulation ist die Standardanwendung der reduction-Anweisung. Es geht darum alle Elemente einer Liste auf irgendeine Weise zu kombinieren. Das Addieren aller Zahlen einer Liste ist ein Beispiel welches mit einer solchen Anweisung gelöst werden kann.

int accumulate(int arr[], int size){
	int sum = 0;
	#pragma omp parallel for reduction(+:sum)
	for(int i=0; i<size; ++i)
		sum += arr;
	return sum;
}

Leider handelt es sich jedoch nicht immer um einen Buildinoperator. Sei [i]Matrix* eine Klasse welche eine 2-dimensionale Matrix beschreibt und für welche der *= Operator entsprechend überladen wurden. Der Defaultkonstruktor initialisiert die Matrix auf die Einheitsmatrix.

class Matrix{
private:
	int 
		a,b,
		c,d
	;
public:
	Matrix():a(1), b(0), c(0), d(1){}
	Matrix&operator*=(const Matrix&right){
		Matrix left = *this;
		a = left.a*right.a + left.b*right.c;
		b = left.a*right.b + left.b*right.d;
		c = left.c*right.a + left.d*right.c;
		d = left.c*right.b + left.d*right.d;
		return *this;
	}
};

Matrix multiply(Matrix arr[], int size){
	const unsigned thread_count = omp_get_max_threads();
	Matrix*sub_prod = new Matrix[thread_count]; 
	#pragma omp parallel num_threads(thread_count)
	{
		const unsigned id = omp_get_thread_num();
		const int begin = size*id/thread_count;
		const int end = size*(id+1)/thread_count;

		for(int j=begin; j<end; ++j)
			sub_prod[id] *= arr[j];
	}
	Matrix prod;
	for(int j=0; j<thread_count; ++j)
		prod *= sub_prod[j];
	delete[]sub_prod;
	return prod;
}

Die Idee ist die Assoziativität der Multiplikation auszunutzen um das Produkt in mehrere Teilprodukte zu zerlegen und diese zuerst von jedem Thread einzeln auswerten zu lassen.

Zuerst wird bestimmt wie viele Threads überhaupt zur Verfügung stehen. Anschießend wird ein Array angelegt welches einen Eintrag für jeden Thread besitzt. In diesem Array werden die Zwischenresultate abgelegt. Die Matrizenmultiplikation ist nicht kommutativ, daher verwende ich zur Berechnung der Teilprodukte keine OpenMP-for-Schleife. Diese garantiert nämliche keine bestimmte Reihenfolge. Eine ordered OpenMP-for-Schleife ist zu steif da nicht alle Multiplikation sequentiell ablaufen müssen. Nachdem alle Teilprodukte berechnet wurden werden sie multipliziert und das Ergebnis wird zurückgegeben.

Ich habe keinen STL-vector verwendet da dieser in der Regel nicht dafür ausgelegt ist, mit mehreren Threads klar zu kommen. In diesem speziellen Fall ist es jedoch höchst unwahrscheinlich, dass dies ein Problem darstellen würde. Man könnte natürlich einen Smartpointer verwenden um die Speicherfreigabe zu garantieren. Allerdings ist das im Fall einer Ausnahme das kleinere Problem.

Die einzigen Operationen die eventuell eine Ausnahme werfen könnten, sind die Speicherallokierung selbst und der *= Operator. Wenn die Speicherallokierung fehlschlägt, geht kein Speicher verloren. Der *= Operator unserer Matrixklasse wirft natürlich keine Ausnahmen allerdings ist dies bei einer allgemeinen Matrixklasse nicht mehr selbstverständlich. In diesem Fall muss man allerdings sicher stellen, dass die Ausnahme sicher aus dem parallel-Block hinaus propagiert wird. Dies ist ein wesentlich größeres Problem als das kleine Speicherleck.

7.1 STL

Im Forum kommt es immer wieder zu Diskussionen zum Thema Parallelisierung und wie man die einsetzt ohne Programme von Grund auf neu zu bauen. Oft kommt dann die Idee auf die Algorithmen der STL zu parallelisieren und den Benutzer-Code unverändert zu lassen. Dies ist allerdings schwer um zu setzen. Im folgenden werde ich zeigen wie man std::accumulate parallelisieren kann und welche Probleme es dabei gibt.

Es sei folgende recht einfache Funktion gegeben:

template<class Iter, class T>
T accumulate(Iter begin, Iter end, T sum){
	for(; begin!=end; ++begin)
		sum += *begin;
	return sum;
}

Unser Ziel ist es diese Funktion durch Benutzung von OpenMP zu beschleunigen. Auf den ersten Blick erscheint dies wie eine sehr einfache Aufgabe, dies ist sie aber nicht!

Zuerst müssen wir sicherstellen, dass die zu Grunde liegenden +-Operation parallel ausgeführt werden kann. Speichert der Operator zum Beispiel Zwischenwerte in einem statischen Array, dann können wir eine Parallelisierung ohne die Operation selbst zu verändern bereits vergessen. Es ist auch möglich, dass Werte vorberechnet werden müssen und der Operator auf diese zugreift, das heißt es muss etwas initialisiert werden. Zum Beispiel eine Primzahlliste oder das Pascalsche Dreieck zur Bestimmung großer Binomialkoeffizienten sind gute Kandidaten. Dies muss vor dem ersten Aufruf passiert sein. Es darf nicht erst bis zur ersten Verwendung gewartet werden ehe man Initialisierungsarbeit erledigt.

Die +-Operation darf keine Ausnahme werfen, da C++ leider keine Möglichkeit bietet diese über Threadgrenzen hinaus zu propagieren.

Des Weiteren benötigt man die Assoziativität der +-Operation. Die Assoziativität ist notwendig damit die Summe in Teilsummen aufgespalten werden kann. Die Kommutativität ist notwendig, da nicht klar ist in welcher Reihe die Threads fertig werden. Dies erscheint offensichtlich, allerdings ist das eine zusätzliche Einschränkung gegenüber der sequentiellen Version. Oft ist sie auch nicht gegeben. Zum Beispiel ist die +-Operation von Gleitkommazahlen nicht assoziativ da die Genauigkeit von der Auswertungsreihenfolge abhängt. Normalerweise sagt man bei Gleitkommazahlen jedoch, dass die Genauigkeit nicht wichtig ist und nimmt einfach an, dass die Assoziativität gegeben ist.

Ohne die Kommutativität kann man auskommen, wenn es sein muss wie das vorherige Beispiel zeigt. Wenn sie allerdings zur Verfügung steht wird die Implementierung einfacher. Ich gehe im Folgenden davon aus, dass sie gegeben ist.

C++-Templates bieten leider praktisch keine Möglichkeiten diese Bedingungen sicher zu stellen. Die Dokumentation der Funktion ist also die einzige Stelle wo man sie angeben kann. Es bleibt nur zu hoffen, dass der Benutzer der Funktion die Dokumentation auch gründlich ließt und versteht...

Ist sicher gestellt, dass die +-Operation parallel ausgeführt werden kann müssen wir zusehen, dass wir wahlfreien Zugriff auf die Liste von Elementen haben. Dies ist nötig da es einem OpenMP-for frei steht in welcher Reihenfolge er die Schleifenduchläufe ausführen will. Ist dies nicht gegeben, kann man entweder die sequentielle Version benutzen oder die Elemente zu erst kopieren. Diese Kopie muss sequentiell erfolgen, da es keinen wahlfreien Zugriff gibt. Was besser ist hängt wieder von der +-Operation ab.

Bei Ganzzahlen ist die sequentielle Version deutlich überlegen, da die Addition wahrscheinlich so schnell ist als das Kopieren einer Zahl. Das einfache Ausrechnen ist also in etwa genauso schnell wie das Vorbereiten der parallelen Berechnung. Bei der Matrizenmultiplikation ist es allerdings besser zu kopieren. Werden Matrizen unterschiedlicher Dimensionen multipliziert so lohnt es sich unter Umständen viel mehr, zuerst dynamisch die optimale Auswertungsreihenfolge zu bestimmen als einfach nur zu parallelisieren. Ich werde im Folgenden davon ausgehen, dass es sich lohnt zu kopieren.

Der Code sieht bis hierhin in etwa folgendermaßen aus:

template<class Iter, class T>
T accumulate(Iter begin, Iter end, T init){
	return accumulate_check_category(begin, end, init, std::iterator_traits<Iter>::iterator_category());
}

template<class Iter, class T>
T accumulate_check_category(Iter begin, Iter end, T init, std::random_access_iterator_tag){
	return accumulate_parallel(begin, end-begin);
}

template<class Iter, class T>
T accumulate_check_category(Iter begin, Iter end, T init, std::input_iterator_tag){
	typename std::iterator_traits<Iter>::value_type T;
	std::vector<T>buffer(begin, end);
	return accumulate_parallel(buffer.begin(), buffer.size());
}

template<class Iter, class T>
T accumulate_parallel(Iter begin, int len, T init){
	// ... 
}

Die eigentliche Parallelisierung ist ein typischer Reduktionsfall. Die entsprechende OpenMP-Anweisung kann aber nur im Fall eines buildin Typen verwendet werden. Da durch die Parallelisierung Kosten entstehen lohnt es sich nicht unter einer Gewissen Grenze zu parallelisieren. Wie viele Element es geben muss, damit sich das Ganze lohnt hängt von der +-Operation, vom Prozessor und dem verwendeten Compiler ab.

Die OpenMP-Anweisung dürft bei einigen Compilern schneller sein als ihr Nachbau. Beim Nachbau muss man zwischen Nachbau mit Kommutativität und ohne Kommutativität unterscheiden. Bei der OpenMP-Anweisung braucht man dies nicht, da es kein nicht kommutativer buildin +-Operator gibt. Dies macht also mindestens 3 Konstanten, welche man praktisch nur empirisch bestimmen kann, die man in den Code einbauen muss. Da ich die Kommutativität angenommen habe sind es bei mir nur 2. Ich nehme einfach mal 50000 für beide an.

Ein weiteres Problem das sich stellt ist die Frage ob die Laufzeit des Operator+ von seinen Operanten abhängt oder nicht. Ob die for-Anweisung ein schedule(dynamic) benötigt oder nicht hängt davon ab. Will man alle Fälle abdecken so muss man jetzt schon 6 Konstanten bestimmen. Ich gehe mal davon aus, dass es in etwa immer gleich schnell ist. Ein Gegenbeispiel wäre die Multiplikation von Matrizen unterschiedlicher Dimensionen.

template<class Iter, class T>
T accumulate_parallel(Iter begin, int len, T init){
	accumulate_parallel_reduction(begin, len, init, boost::is_arithmetic<typename std::iterator_traits<T>::value_type>::type());
}

template<class Iter, class T>
T accumulate_parallel_reduction(Iter begin, int len, T sum, boost::true_type){
	if(len < 50000){
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}else{
		#pragma omp parallel reduction(+:ret) for
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}
	return sum;
}

template<class Iter, class T>
T accumulate_parallel_reduction(Iter begin, int len, T sum, boost::false_type){
	if(len < 50000){
		for(int i=0; i<len; ++i)
			sum += *(begin+i);
	}else{
		#pragma omp parallel 
		{
			T sub_sum = T();
			for(int i=0; i<len; ++i)
				sub_sum += *(begin+i);
			#pragma omp critical
			{
				#pragma omp flush(sum)
				sum += sub_sum;
				#pragma omp flush(sum)
			}
		}
	}
	return sum;
}

Ich habe auf meinem 4-Kern-Prozessor und dem GCC-4.2 ein paar Tests gemacht um zu bestimmen wie groß bei mir diese Konstanten liegen. Unter etwa 50000 int-Elementen ist die sequentielle Version bei mir im Vorteil. Bei mehr Elementen ist die parallelisierte Version besser. Erstaunlich ist, dass die Reduktions-Answeisung in etwa genau so schnell ist wie ihr Nachbau. Ich gehe mal davon aus, dass die GCC-Entwickler den OpenMP Code noch nicht wirklich optimiert haben und in dieser Version vor allem auf korrektes Verhalten Wert gelegt haben.

8 Gibt es noch mehr?

In diesem Artikel habe ich nicht alle OpenMP-Anweisungen vorgestellt. Ich habe mich legendlich auf die mir am wichtigsten erscheinenden beschränkt. Wer OpenMP wirklich produktiv einsetzen will, dem empfehle ich die Spezifikation durch zu lesen. Für eine Spezifikation ist diese erstaunlich leicht verständlich und es sind auch viele Beispielanwendungen darin enthalten.

9 Quellen

• http://bisqwit.iki.fi/story/howto/openmp/
• http://www.openmp.org/
• http://algo2.iti.uni-karlsruhe.de/singler/mcstl/