© alphaspirit, 123RF

Besonders im Bereich der Natur- und Ingenieurwissenschaften ist Parallelverarbeitung heute unverzichtbar. Aber auch der normale Desktopanwender kann mit mindestens vier Rechenkernen von höherer Performance durch parallele Ausführung profitieren.

Open MP, Threading oder MPI – viele Programmiersprachen bieten Konzepte für die Parallelverarbeitung an. Diese Spracherweiterungen sind meist recht einfach, aber ihre effiziente Anwendung ist schwierig, weil sich nicht alle Algorithmen dafür einfach umschreiben lassen. Oft sind komplett neue Ansätze notwendig. Auch High-Level-Sprachen wie Python oder R offerieren gewisse parallele Erweiterungen. Allerdings sind diese erstens nicht von Anfang an im Konzept enthalten (intrinsisch) und zweitens sind Sprachen wie Python und R oft extrem langsam, wenn es um numerische Berechnungen geht.

Es gibt also eine Marktlücke für eine Programmiersprache, die intrinsisch parallel, einfach zu handhaben und gleichzeitig noch sehr performant ist. Julia versucht genau diese Lücke zu schließen. Tabelle 1 zeigt einige auf der Julia-Website veröffentlichte Micro-Benchmarks für verschiedene Programmiersprachen und Probleme. Die Werte sind so normalisiert, dass C den Wert 1 hat und kleinere Werte bessere Performance bedeuten. Für einige der Micro-Benchmarks erreicht Julia fast C-Performance und ist viel schneller als Sprachen wie R und Python.

Die Performance lässt sich auch an noch einfacheren Beispielen demonstrieren. Dazu soll sowohl Python als auch Julia viele Quadratwurzeln in einer Schleife berechnen. Die Beispiele schreiben die Schleife explizit aus:

import math
import time
tstart = time.time()
for i in xrange(100000000):
 math.sqrt(i)
tstop = time.time()
print tstop-tstart

Auf dem Testsystem benötigt diese Schleife zirka 10 Sekunden. Exakt die gleichen Operationen führt der folgende Julia-Code aus:

tstart = time()
for i = 1:100000000
 sqrt(i);
end
tstop = time()
println(tstop-tstart)

Dieser Code braucht weniger als 0,1 Sekunden, ist also um mehr als den Faktor 100 schneller als der äquivalente Python-Code. Hier liegt eine der Hauptstärken von Julia: Multiple Dispatch bei Funktionsaufrufen führt zu extrem effizienten Code, der praktisch jeder High-Level-Sprache in puncto Geschwindigkeit überlegen ist. Schneller Code in Julia lässt sich ohne Tricks wie Vektorisierung oder Auslagerung in C-Erweiterungen erreichen. Dagegen sind solche Tricks häufig nötig, wenn Python- oder R-Code zu beschleunigen ist.

Hello Julia!

Die aktuelle Julia-Version ist über Github zu beziehen:

git clone https://github.com/JuliaLang/julia.git

Auch die meisten Linux-Distributionen enthalten fertige Julia-Pakete, die ebenfalls verwendbar sind. Julia ist aber eine sehr junge Sprache, die stetig weiterentwickelt wird. Es empfiehlt sich daher, immer mal einen Blick auf die aktuelle Entwickler-Release zu werfen.

Wer sich direkt für die Github-Version entscheidet, muss noch den Sourcecode via »make« übersetzen. Das kann durchaus einige Minuten dauern, da der Quelltext recht umfangreich ist. Auf den meisten Systemen sollte alles problemlos laufen, allerdings erklärt die Datei »README.md« einige Fälle, in denen manuell nachzujustieren ist. Nach dem Übersetzen findet sich die Binary-Datei »julia« im Hauptverzeichnis. Wer sie startet, der sieht den in Abbildung 1 gezeigten Julia-Prompt.

Abbildung 1: Der Begrüßungsbildschirm von Julia direkt nach dem Start.

Hier beginnt der Reigen mit einem Befehl. Guter alter Tradition folgend kann der Programmierer ein Hello-World-Beispiel probieren:

julia> for i=1:10
 println("Hello World ",i)
 end

ergibt

Hello World 1 Hello World 2 Hello World 3 ...

Die Syntax von Julia ist stark an die von Matlab angelehnt. So beginnt beispielsweise die Indizierung von Arrays bei 1 und nicht bei 0. Auch die Syntax von For-Schleifen ist exakt gleich der von Matlab. Der Fokus soll im Folgenden aber auf den parallelen Fähigkeit von Julia liegen. Die grundlegenden Sprachkonstrukte erklärt die offiziellen Julia-Dokumentation [1] ausführlich.

Parallel mit Julia

Es gibt mehrere Arten paralleler Programmierung. Vielen Anwendungsprogrammierern dürfte beispielsweise Threading geläufig sein. Damit lassen sich sehr einfach Schleifen parallelisieren. Dazu bricht der Programmierer die Schleife auf und ordnet jedem Prozessor des Systems einen Teil der Schleife zu. Da alle Prozessoren parallel rechnen, wird die gesamte Schleife schneller abgearbeitet.

Besonders im High Performance Computing ist jedoch eine andere Art der Parallelisierung weit verbreitet: Distributed-Memory-Parallelisierung. Hier wird nicht nur die Rechenarbeit verteilt, sondern auch die Daten. Häufig kommt dabei Message-Passing zum Einsatz, wie es beispielsweise im MPI-Standard implementiert ist.

Julia verwendet noch einen anderen Ansatz, nämlich eine Master-Worker-Architektur oder One-Sided-Message-Passing. Ein Prozess (der Master) gibt hierbei Anweisungen an andere Prozessoren (Sklaven). Der Programmierer muss also nur explizit einen Prozess verwalten. Die gesamte Julia-Parallelisierung ist auf zwei Hauptkonstrukten aufgebaut: Remote References und Remote Calls.

Eine Remote Referenz verweist auf ein Objekt, das mit einem anderen Prozess verbunden ist. Beispielsweise lässt sich das Ergebnis einer Rechnung, die ein anderer Prozess erbracht hat, lokal verwenden, indem man eine solche Referenz benutzt. Genauso ist ein Remote Call ein Aufruf, der eine Funktion von einem anderen Prozess ausführen lässt. So darf der Programmierer beliebige Julia-Funktionen via Remote Call von beliebigen Prozessen ausführen lassen.

Das ermöglicht auf einfache Weise verteiltes Rechnen mit Julia. Alle anderen parallelen Sprachkonstrukte in Julia verwenden diese beiden Konzepte, um parallele Aufgaben zu bewältigen.

Um Julia parallel zu betreiben, muss der Anwender dem Programm beim Start mitteilen, wie viele Prozesse zu starten sind. Hat die Maschine zum Beispiel vier Prozessor-Cores, so könnte er vier Prozesse via

./julia -p 4

starten. Es wäre auch möglich, mehr Prozesse zu starten, allerdings wird dies bei nur vier physikalisch vorhandenen Cores keinen Performancegewinn bringen.

Prozesse lassen sich mit der »addprocs()« -Funktion auch on the Fly hinzufügen. Jeder Prozess hat eine ID, wobei dem ersten Masterprozess die ID 1 zugeordnet ist. Alle Worker-Prozesse haben IDs größer als 1. Die »myid()« -Funktion gibt die ID eines Prozesses zurück.

Julia kann auch parallel auf ganzen Clustern von Rechnern laufen. Dazu muss eine passwortlose SSH-Umgebung vorhanden sein. Mit der »machinefile« -Option fügt Julia dann die einzelnen Clusterrechner hinzu.

Wer eine Aufgabe von einem bestimmten Prozess ausführen lassen will, verwendet dazu die Low-Level-Funktion »remotecall()« :

julia> result = remotecall(3, +, 2, 2)
RemoteRef(3,1,6)

Der Prozess mit der ID 3 soll 2 + 2 berechnen, ein triviales Beispiel. Als Antwort kommt allerdings nicht direkt das Ergebnis, sondern eine Remote Reference (»RemoteRef« ), denn das Ergebnis ist nicht lokal verfügbar, da es der Prozess 3 berechnete. Um das Ergebnis zu erhalten, muss der es erst von dem entfernten Prozess via »fetch()« abholen:

julia> fetch(result)
4

Wie bereits erwähnt, ist »remotecall()« eine Low-Level-Funktion. Ist es egal, welcher Prozess die Addition ausführt, so ist Julias »@spawn« -Makro verwendbar, um Aufgaben an andere Prozesse zu senden. Das »@spawnat« -Makro erlaubt es zusätzlich noch, den Prozess genau zu spezifizieren. So erzeugt

julia> @spawnat 2 rand(10,10)

eine Zufallsmatrix durch Prozess 2, während

julia> @spawn rand(10,10)

Julia die Wahl des Prozesses überlässt. Je nach Algorithmus und Problem kann der Programmierer also entscheiden, ob er die Prozessverteilung selbst in die Hand nimmt oder dies komplett Julia machen lässt. In den meisten Fällen ist es am einfachsten, diese Aufgabe an Julia zu delegieren. Auch »@spawn« liefert lediglich eine Remote Reference zurück, die dann wieder einen »fetch()« -Aufruf zum Abholen des Ergebnisses benötigt:

julia> result = @spawn 2+2
RemoteRef(2,1,6)julia> fetch(result)
4

Ein Aufruf von »remotecall()« oder von »@spawn« garantiert nicht, dass die unmittelbar zurückgelieferte Remote Reference auch sofort schon das Ergebnis enthält. Erst »fetch()« stellt dies sicher. Beispielsweise könnte »@spawn« eine komplexe und langwierige Rechnung mit Hilfe eines anderen Prozesses starten. Trotzdem wird der »@spawn« -Aufruf zunächst nicht blockieren.

Allerdings kann»fetch()« erst dann das Ergebnis liefern, wenn die Berechnung auch abgeschlossen ist. Das heißt, es ist »fetch()« , das dann eine Zeit lang blockiert und der lokale Prozess wartet so lange, bis das Ergebnis verfügbar ist.

Verteilte Objekte

Julia bietet auch die Möglichkeit, verteilte Objekte direkt zu generieren. Das ist insbesondere beim Arbeiten mit größeren Matrizen hilfreich. Die Sprache stellt hierzu so genannte verteilte Arrays zur Verfügung, die den Typ »DArray« haben. Um diese Arrays nutzen zu können, muss der Anwender zunächst ein entsprechendes Paket installierten. Das geht in Julia sehr einfach mit Hilfe der internen Paketverwaltung (»Pkg« ):

julia> Pkg.checkout("DistributedArrays")

Damit wird der »DistributedArrays« -Code heruntergeladen und installiert. Eine Übersicht aller installierter Pakete ist via »Pkg.status()« möglich. Sollten Updates von Paketen verfügbar sein, so lassen die sich einfach via »Pkg.Update()« einspielen. Um verteilte Arrays nutzen zu können, muss noch sichergestellt werden, dass alle Prozesse das Paket laden:

julia> @everywhere using DistributedArrays

Das »@everywhere« -Makro lädt das Paket via »using« und sorgt dafür, dass es allen Prozessen zur Verfügung steht. Grundsätzlich wird eine Anweisung hinter »@everywhere« von jedem Prozess aufgerufen.

Es gibt zahlreiche Möglichkeiten, verteilte Arrays zu erzeugen. Am einfachsten geschieht dies mit folgenden Funktionen:

dzeros(100,100,10)
dones(100,100,10)
drand(100,100,10)
drandn(100,100,10)
dfill(x,100,100,10)

Wie die Namen schon ahnen lassen, stellen die Funktionen verteilte Arrays mit den angegebenen Dimensionen bereit, die beispielsweise mit Nullen gefüllt sind (»dzeros« ). Julia übernimmt hier die gesamte Arbeit: Die einzelnen Matrixelemente sind anschließend auf die verschiedenen Prozesse verteilt und auch initialisiert.

Julia stellt auch einige hilfreiche Funktionen bereit, um mit verteilten Arrays zu arbeiten. Die »localpart()« -Funktion gibt beispielsweise Teil jenes Array wieder, das dem lokalen Prozess zugewiesenen ist. Analog gibt »localindexes()« die Indizes des lokalen Teils des Array zurück. Soll der lokale Prozess ein verteiltes Array komplett bearbeiten, so kann der Programmierer dies mit »convert()« bewerkstelligen. Vice versa wandelt »distribute()« ein lokales Array in ein verteiltes Array um, das heißt, die Funktion verteilt die Daten des Array auf die vorhandenen Prozesse.

Initialisierung

Neben den oben genannten Funktionen zum Anlegen von Arrays lassen sich auch komplizierte Initialisierungsfunktionen verwenden. Dazu dient:

DArray(init, dims[, procs, dist])

Die Init-Funktion muss der Programmierer bereitstellen. Sie bekommt ein Tupel von Index-Bereichen übergeben, für die das Array dann innerhalb der Init-Funktion initialisiert wird – die Init-Funktion initialisiert also lokale Bereiche des verteilten Array. Die »dims« -Variable gibt wie üblich die Dimensionen des Array an. Optional lässt sich via »procs« noch bestimmen, über welche Prozesse das Array zu verteilen ist.

Außerdem kann der Programmierer mit »dist« noch festlegen, wie das Array auf die verschiedenen Prozesse zu verteilen ist. Das ist aber im Allgemeinen gar nicht nötig, da Julia das automatisch regelt. Das Beispiel in Listing 1 erzeugt eine verteilte Zufallsmatrix und gibt dann jenen Teil zurück, für den Prozess 2 verantwortlich ist.

01 julia> @everywhere using DistributedArrays
02
03 julia> d=drand(10,10);
04
05 julia> r=@spawnat 2 begin
06  localpart(d)
07  end
08
09 RemoteRef(2,1,23)
10
11 julia> fetch(r)
12 5x5 Array{Float64,2}:
13  0.932366 0.181326 0.0261475 0.211363 0.308226
14  0.330008 0.924271 0.543386 0.895825 0.617452
15  0.51471 0.801718 0.786854 0.174585 0.413264
16  0.840219 0.750775 0.126154 0.853618 0.899762
17  0.866529 0.804654 0.19808 0.49411 0.951355

Das Praktische an diesen verteilten Arrays ist, dass Julia die gesamte Kommunikation vor dem Benutzer versteckt. Will der beispielsweise die Summe aller Elemente der »d« -Matrix berechnen, so ruft er einfach »sum(d)« auf. Das liefert direkt das Ergebnis. Julia summiert dabei mit jedem einzelnen Prozess alle Elemente des Teilarray auf. Anschließend werden dann alle Teilergebnisse zur Gesamtsumme addiert.

Julia bietet auch Konstrukte, um Schleifen einfach zu verteilen. Folgende wenige Zeilen addieren beispielsweise mit »rand()« generierte Zufallszahlen:

julia> r = @parallel (+) for i=1:200000000
 rand()
 end

Das »@parallel« -Makro verteilt die Schleife auf alle zur Verfügung stehenden Prozesse. Der Rückgabewert eines jeden Schleifendurchlaufs ist dabei der letzte Ausdruck in der Schleife, also hier die »rand()« -Anweisung zum Erzeugen einer Zufallszahl. Die Ergebnisse eines jeden Schleifendurchlaufs addiert dann das obige Beispiel. Den Reduktionsoperator (hier die Addition) kann der Programmierer auch weglassen. Dann verteilt Julia die Schleife parallel, ohne die Ergebnisse am Ende zu reduzieren.

Für solche Fälle bietet sich auch die »pmap()« -Funktion an. Sie führt einfach eine Funktion eines bestimmten Objekts verteilt aus, also eine typische Map-Aufgabe. Folgendes Beispiel nutzt die »pmap()« -Anweisung, um den Rang von Matrizen zu berechnen:

julia> M = [rand(1000,1000) for i=1:4];julia> pmap(Rang,M)
4-element Array{Any,1}:
 1000
 1000
 1000
 1000

Die erste Anweisung erzeugt vier 1000-mal-1000-Matrizen. Dann berechnet »pmap()« den Rang jeder Matrix verteilt. Da Julia mit »-p 4« gestartet wurde, ist jeweils ein Worker-Prozess für eine Matrix zuständig. Genauso einfach lassen sich verteilt andere Eigenschaften der Matrizen berechnen. Zum Beispiel berechnet die »det()« -Funktion die Determinante, die »inv()« -Funktion die Inverse einer Matrix.

Monte-Carlo-Berechnungen lassen sich typischerweise sehr einfach parallelisieren, da die einzelnen Rechenschritte meist unabhängig voneinander sind.

Ein erstes Beispiel: Monte-Carlo-Berechnung von Pi

Dies ist auch der Fall bei der Monte-Carlo-basierten Berechnung der Kreiszahl Pi. Dazu generiert der Programmierer sehr viele Zufallszahlen in einem den Einheitskreis (Radius 1) umgebenden Quadrat. Diese Zufallszahlen werden gleichmäßig zwischen -1 und +1 für die x- und y-Koordinate erzeugt.

Das Verhältnis der Fläche des Einheitskreises zu der des Quadrats beträgt nun gerade pi*1*1/(2*2). Das bedeutet, dass genau der Anteil pi/4 der Zufallszahlen im Kreis liegen sollte. Zählt man die nun in einer Schleife, so ergibt sich als Monte-Carlo-Schätzung für Pi:

pi_MC = N_in / N_tot * 4

Dabei ist »N_in« die Zahl der Zufallszahlen im Einheitskreis und »N_tot« die Gesamtzahl der generierten Zufallszahlen in der x-y-Ebene. Je höher »N_tot« , umso näher rückt »pi_MC« an Pi heran. Es ist also wünschenswert, viele Zufallszahlen zu generieren, um einen möglichst genauen Wert für Pi zu erhalten.

Listing 2 zeigt eine Implementierung dieses Verfahrens in Julia. Sie führt die zentrale Schleife mit »@parallel« verteilt aus und summiert das Schleifenergebnis (1 oder 0) am Ende in einem Reduktionsschritt auf. Auf dem Testsystem läuft dieser Code auf einem Prozessor in knapp zwei Minuten:

01 N_tot = 1000000000
02
03 tstart = time()
04
05 N_in = @parallel (+) for n=1:N_tot
06         x = rand() * 2 - 1
07         y = rand() * 2 - 1
08
09         r2 = x*x + y*y
10         if r2 < 1.0
11                 1
12         else
13                 0
14         end
15 end
16
17 tstop = time()
18
19 pi_MC = N_in/N_tot*4.0
20
21 println("time        = ", tstop-tstart, " seconds")
22 println("pi estimate = ", pi_MC)

time = 116.90136814117432 seconds
pi estimate = 3.14162308

Wer zwei Cores nutzt, beschleunigt das Beispiel fast um den Faktor 2:

time = 73.63914084434509 seconds
pi estimate = 3.141607444

Das Verfahren konvergiert recht langsam, aber die ersten beiden Stellen von Pi (3,14) sind immerhin schon korrekt berechnet.

Ein zweites Beispiel: Julia mit Julia

Als ein etwas komplizierteres Beispiel folgt hier die parallele Berechnung eines Julia-Fraktals mit Julia. Das Programm hat ursprünglich Martin Rupp basierend auf dem offiziellen Julia-Mandelbrot-Beispiel geschrieben [2]. Die hier gezeigte Version ist allerdings modifiziert, da sich die Syntax verändert hat.

Abbildung 2: Julia berechnet Julia-Mengen mit verteilten Arrays.

Die mathematische Idee eines Julia-Fraktals ist sehr einfach. Man betrachtet dabei die Folge komplexer Zahlen zn+1 = zn + c für eine beliebige komplexe Konstante c. Jedes c ergibt eine bestimmte Julia-Menge. Die typischen Julia-Bilder entstehen durch Einsetzen eines Zählers. Je mehr Iterationen in der zn-Folge zum Erreichen eines bestimmten Schwellenwertes nötig sind, desto heller erscheint der zugehörige Punkt. Die zweidimensionale komplexe Ebene lässt sich einfach auf die x- und y-Koordinate eines Pixelrasters abbilden, wodurch die bekannten Bilder entstehen.

01 #DistributedArrays, WIDTH, HEIGHT und MAXITER sind global
02 @everywhere using DistributedArrays
03 @everywhere WIDTH=1000
04 @everywhere HEIGHT=500
05 @everywhere MAXITER=100
06
07 #Images ist lokal
08 using Images
09
10
11
12 # Julia-Funktion
13 @everywhere function julia(z, maxiter::Int64)
14  c=-0.8+0.156im
15  for n = 1:maxiter
16  if abs(z) > 2.0
17  return n-1
18  end
19  z = z^2 + c
20  end
21  return maxiter
22 end
23
24
25
26 # Init-Funktion fuer das Anlegen des Array
27 @everywhere function parJuliaInit(I)
28  # lokaler Patch
29  d=(size(I[1], 1), size(I[2], 1))
30  m = Array(Int, d)
31
32  xmin=I[2][1]
33  ymin=I[1][1]
34
35  for x=I[2], y=I[1]
36  c = complex((x-WIDTH/2)/(HEIGHT/2), (y-HEIGHT/2)/(HEIGHT/2))
37  m[y-ymin+1, x-xmin+1] = julia(c, MAXITER)
38  end
39  return m
40 end
41
42
43
44 # verteiltes Array anlegen
45 Dm = DArray(parJuliaInit, (HEIGHT, WIDTH))
46
47 # verteiltes Array auf den lokalen Prozess holen
48 m = convert(Array, Dm)
49
50 #Bild als PNG speichern
51 imwrite(grayim(transpose(m)/(1.0*MAXITER)),"juliaset.png")

Der in Listing 3 abgedruckte Code macht genau das. Julia kann problemlos mit komplexen Zahlen umgehen, weshalb die Implementierung besonders einfach ist. Kompliziert wird der Code durch die Parallelisierung, weil verteilte Arrays zum Einsatz kommen. Je nach Prozess-ID wird so jeweils ein Teil des Julia-Bildes von einem anderen Prozess berechnet. Ein Vorteil dieser Parallelisierung ist, dass die Berechnung eines bestimmten Pixels unabhängig von dessen Umgebung ist. Komplizierter wäre die Situation dann, wenn die Eigenschaft eines Pixels auch von der Umgebung des Pixels abhinge. So lässt sich mit

~/julia/julia -p 4 juliaset.jl

der Code einfach ausführen (siehe Abbildung 2).

Fazit

Obwohl Julia noch eine junge Sprache ist, kann man mit ihr schon viel produktiven Code schreiben. Gerade für Anwendungen im Bereich numerischer Berechnungen ist Julia derzeit praktisch unschlagbar und wird ältere Toolsets wie Python, R oder Matlab ablösen.