Eine Parallelisierung der schnellen Fourier

Eine Parallelisierung der schnellen FourierTransformation auf SO(3)
A Parallelization of the Fast Fourier Transform on SO(3)
Bachelorarbeit
Im Rahmen des Studiengangs
Informatik
Vorgelegt von
Denis-Michael Lux
Ausgegeben und betreut von
Prof. Dr. Jürgen Prestin
der Universität zu Lübeck
mit Unterstützung von
Christian Wülker, M.Sc.
Lübeck, den 06. Oktober 2015
Selbstständigkeitserklärung
Der Verfasser versichert an Eides statt, dass er die vorliegende Arbeit selbständig, ohne fremde Hilfe und ohne Benutzung anderer als der angegebenen Hilfsmittel angefertigt hat.
Die aus fremden Quellen (einschließlich elektronischer Quellen) direkt oder indirekt übernommenen Gedanken sind ausnahmslos als solche kenntlich gemacht. Die Arbeit ist in gleicher oder ähnlicher Form oder auszugsweise im Rahmen einer anderen Prüfung noch nicht
vorgelegt worden.
Ort/Datum
Unterschrift des Verfassers
Inhaltsverzeichnis
Einleitung
1
1 Mathematische Theorie
3
1.1 Rotationsgruppe SO(3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2 Wigner-D-Funktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.3 Diskrete Fourier-Transformation auf SO(3) (DSOFT) und diskrete inverse
Transformation (iDSOFT) . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
1.4 Schnelle Fourier-Transformation auf SO(3) (SOFFT) und schnelle inverse
Transformation (iSOFFT) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 Implementierung und Parallelisierung
10
13
2.1 Sequentielle SOFFT und iSOFFT . . . . . . . . . . . . . . . . . . . . . . . .
13
2.2 Open Multi-Processing (OpenMP) . . . . . . . . . . . . . . . . . . . . . . .
16
2.3 Entwurf paralleler Algorithmen mittels Forster-Methodik . . . . . . . . . .
17
2.4 Parallelisierung der SOFFT/iSOFFT mittels Forster-Methodik und OpenMP
19
2.5 Analyse der Parallelisierung . . . . . . . . . . . . . . . . . . . . . . . . . .
25
3 Numerische Beispiele
33
3.1 Benchmark zu Speedup, Eﬃzienz und Rechengenauigkeit . . . . . . . . . .
33
3.2 Ergebnisse der Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
4 Diskussion und Schlussfolgerung
41
Abkürzungsverzeichnis
45
Algorithmenverzeichnis
47
Abbildungsverzeichnis
49
Literaturverzeichnis
51
Zusammenfassung
In dieser Arbeit wird eine Parallelisierung der schnellen Fourier-Transformation auf
der Rotationsgruppe SO(3) (SOFFT) und ihrer Inversen (iSOFFT) entworfen, implementiert und in der Theorie und Praxis analysiert. Zu diesem Zweck werden zunächst
die mathematischen Grundlagen dieser Algorithmen vorgestellt. Anschließend wird eine sequentielle Implementierung der SOFFT und iSOFFT präsentiert und besprochen.
Für den Entwurf und die Umsetzung der Parallelisierung dieser sequentiellen Implementierung wird die Forster-Methodik und der OpenMP-Standard verwendet. Die Arbeitsoptimalität der entstandenen Parallelisierung wird daraufhin nachgewiesen. Abschließend wird die Praxistauglichkeit der parallelisierten SOFFT und iSOFFT anhand
dreier Benchmarks demonstriert.
Abstract
In this work, a parallelization of the fast Fourier transform on the rotation group SO(3)
(SOFFT) and its inverse (iSOFFT) is designed, implemented, and analyzed in theory
and practice. For this purpose, the underlying mathematical theory of these algorithms is introduced first. Subsequently, a sequential implementation of the SOFFT and
iSOFFT is presented and discussed. For the design and implementation of the parallelization of the sequential implementation, the Forster methodology and the OpenMP
standard are used. The work optimality of the resulting parallelization is then proven.
Finally the applicability of the parallelized SOFFT and iSOFFT in practice is demonstrated by means of three benchmarks.
1
Einleitung
In vielen wissenschaftlichen Anwendungen und Problemen spielen Rotationen im Anschauungsraum R3 eine wichtige Rolle. Bei der Technik des sogenannten fast rotational matchings
[Kovacs und Wriggers, 2002] wird beispielsweise versucht, diejenige Rotation zu bestimmen,
unter der ein gegebenes Objekt möglichst gut mit einem weiteren Objekt übereinstimmt. Die
Bedeutung des Begriﬀs der Übereinstimmung zweier Objekte ergibt sich dabei in der Regel
aus dem speziellen Kontext. Betrachtet man zum Beispiel gegeneinander rotierte dreidimensionale Bilddatensätze zweier medizinischer Bildgebungsmodalitäten, stellt sich oft die Frage, wie ähnlich sich diese Bilder sind. Bei der Bildregistrierung kann versucht werden, diese
Ähnlichkeit durch Rotationen eines der beiden Bilder zu vergrößern [Modersitzki, 2009,
Abschnitte 1.1 und 4.4]. Ein anderes Beispiel findet sich im Bereich der Molekularbiologie.
Hier ist es oft das Ziel, den Überlapp der Elektronendichtefunktionen zweier gegebener Makromoleküle zu maximieren. Die Bestimmung des maximalen Überlapps kann dann dazu
genutzt werden, die Ähnlichkeit der beiden Makromoleküle zu beurteilen. Dadurch ist es
möglich, Rückschlüsse auf Übereinstimmungen in den Eigenschaften und Funktionen dieser Moleküle zu ziehen. Dieses Vorgehen bildet den ersten Schritt beim sogenannten virtual
drug screening [Ritchie, 2011, Abschnitt 1.1.7].
Die von Kostelec und Rockmore [2008] vorgestellte schnelle Fourier-Transformation auf der
Rotationsgruppe SO(3) (SOFFT) und ihre Inverse (iSOFFT) bilden heute einen festen Bestandteil in vielen auf dem fast rotational matching basierenden Verfahren (siehe zum Beispiel [Bülow und Birk, 2011a,b, Slabaugh et al., 2008]). Dies ist auch darauf zurückzuführen,
dass für viele Anwendungsbereiche die Praxistauglichkeit der SOFFT und iSOFFT in [Kostelec und Rockmore, 2008] bereits nachgewiesen wurde. Bis heute wurde eine Parallelisierung dieser Algorithmen allerdings noch nicht vorgenommen, beziehungsweise noch nicht
vorgestellt. Eine solche Parallelisierung der SOFFT und iSOFFT zu entwerfen, umzusetzen
und die Praxistauglichkeit der daraus resultierenden Implementierung nachzuweisen, ist das
Ziel dieser Arbeit.
Der Inhalt ist wie folgt gegliedert. Zu Beginn werden in Kapitel 1 die mathematischen Grundlagen der Rotationsgruppe und der darauf definierten Fourier-Transformation vorgestellt.
Daraus wird dann die SOFFT und die iSOFFT abgeleitet. Im Anschluss daran wird in Kapitel 2 zunächst eine sequentielle Implementierung dieser beiden Algorithmen präsentiert.
Nach einer kurzen Einführung in die C-, C++- und Fortran-Spracherweiterung OpenMP
und die sogenannte Forster-Methodik wird eine Parallelisierung der sequentiellen Implementierung entworfen und umgesetzt. In dem darauﬀolgenden Theorie-Abschnitt wird die Arbeitsoptimalität dieser Parallelisierung nachgewiesen. In Kapitel 3 werden dann drei Benchmarks vorgestellt, anhand derer im abschließenden Kapitel 4 der praktische Nutzen der par-
2
Einleitung
allelisierten Implementierung gezeigt wird. Außerdem werden Ansatzpunkte für mögliche
Weiterentwicklungen genannt.
* * * * *
Einen ganz besonderen Dank möchte ich meinem Betreuer Prof. Dr. Jürgen Prestin entgegenbringen, für seine Mühe, sich diese Arbeit durchgelesen und mit vielen hilfreichen Kommentaren versehen zu haben. Ein weiterer besonderer Dank gebührt Christian Wülker für
die vielen Kommentare und oft (sehr) langen Arbeitstage des Korrekturlesens der einzelnen
Kapitel, sowie der zahlreichen Gesprächen über das entsprechend Thema. Ohne die vielen
Kommentare meiner Betreuer würde die Arbeit heute anders aussehen.
Auch meiner Familie sei für die Unterstützung in meinem Studium gedankt.
Des Weiteren danke ich meinen Kommilitonen für zahlreiche anregende und aufschlussreiche Diskussionen. In diesem Zusammenhang möchte ich Thomas Buddenkotte für die vielen
Gespräche zu seiner Mathematica-Implementierung der SOFFT und iSOFFT danken, sowie
Christoph Schulz für den Hinweis auf diese Thematik als mögliche Bachelorarbeit.
3
1
Mathematische Theorie
In diesem Kapitel wird die schnelle Fourier-Transformation und schnelle inverse FourierTransformation auf der Rotationsgruppe eingeführt. Dazu wird in Abschnitt 1.1 zunächst die
Rotationsgruppe definiert und beschrieben. In Abschnitt 1.2 wird dann der Vektorraum der
Funktionen, welche von der Rotationsgruppe in die komplexen Zahlen abbilden, betrachtet.
Von besonderem Interesse sind in diesem Kontext die sogenannten Wigner-D-Funktionen
und die von diesen Funktionen aufgespannten Untervektorräume. Die betrachteten Untervektorräume sind stets endlichdimensional und besitzen ein Skalarprodukt, bezüglich dessen die Wigner-D-Funktionen eine Orthonormalbasis darstellen. Dies erlaubt es, in Abschnitt 1.3 mittels eines geeigneten Abtasttheorems eine diskrete Fourier-Transformation
und die dazugehörige diskrete inverse Fourier-Transformation auf diesen Untervektorräumen zu definieren. In Abschnitt 1.4 wird aus diesen Transformationen abschließend die entsprechende schnelle Fourier-Transformation und schnelle inverse Fourier-Transformation
abgeleitet.
1.1 Rotationsgruppe SO(3)
Definition 1.1 ([Fischer, 2008, S. 4]). Eine Menge G zusammen mit einer Verknüpfung
: G × G → G, (a, b) → a b,
bezeichnet als (G, ), heißt Gruppe, wenn folgendes gilt:
• Die Verknüpfung is assoziativ, d.h. (a b) c = a (b c) für alle a, b, c ∈ G.
• Es existiert genau ein eindeutig bestimmtes 1 ∈ G mit 1 a = a 1 = a für alle
a ∈ G (1 heißt neutrales Element von G).
• Zu jedem a ∈ G gibt es ein eindeutig bestimmtes a−1 ∈ G mit a−1 a = a a−1 = 1
(a−1 heißt Inverses von a).
Anmerkung. Erfüllt eine Gruppe G zu den bereits genannten Bedingungen zusätzlich die
Bedingung der Kommutativität, d.h.
ab=ba
für alle a, b ∈ G, so wird G auch abelsch genannt.
4
Mathematische Theorie
Definition 1.2 ([Fischer, 2008, S. 13]). Die orthogonale Gruppe im dreidimensionalen reellen
Raum O(3) ist definiert als die Menge aller orthogonalen reellen 3×3-Matrizen mit Determinante +1 oder -1, d.h.
O(3) := R ∈ R3×3 | RT R = 1 ,
wobei 1 die 3×3-Identitätsmatrix und die Verknüpfung die Matrix-Multiplikation bezeichnet. Das neutrale Element von O(3) ist die Identitätsmatrix 1 und das inverse Element für
ein beliebiges Element R ∈ O(3) ist R−1 = RT .
Definition 1.3 ([Fischer, 2008, S. 23]). Die spezielle orthogonale Gruppe im dreidimensionalen reellen Raum oder auch Rotationsgruppe SO(3) ist definiert als die Menge aller reellen
3×3-Matrizen mit Determinante +1, d.h.
SO(3) := R ∈ R3×3 | RT R = 1, det(R) = 1 ⊂ O(3).
Die Verknüpfung beschreibt die Matrix-Multiplikation, das neutrale Element von SO(3) ist
die 3×3-Identitätsmatrix 1 und das inverse Element zu einem beliebigem Element R ∈ SO(3)
ist R−1 = RT .
Anmerkung. Aus Sicht der linearen Algebra beschreibt jedes Element R ∈ SO(3) eine Rotation um den Ursprung des dreidimensionalen reellen Raums R3 . Daher bezeichnet im Folgenden eine Rotation stets ein Element R ∈ SO(3). Darüber hinaus erfüllt die Gruppe SO(3)
nicht die Bedingung der Kommutativität und ist somit keine abelsche Gruppe.
Definition 1.4 ([Edmonds, 1996, S. 53]). Die elementaren Rotationsmatrizen Rx (α), Ry (α)
und Rz (α) sind für einen beliebigen Winkel α ∈ R im Bogenmaß definiert als:


1
0
0



Rx (α) := 
0 cos α − sin α,
0 sin α cos α


Ry (α) := 

cos α
0 sin α


0 
,
− sin α 0 cos α
0
1


cos α − sin α 0



Rz (α) := 
 sin α cos α 0,
0
0
1
wobei Rx (α) eine Rotation um die x-Achse, Ry (α) eine Rotation um die y-Achse und Rz (α)
eine Rotation um die z-Achse repräsentiert.
Wigner-D-Funktionen
5
Satz 1.5 ([Kostelec und Rockmore, 2008, S. 156]). Für jede Rotation R ∈ SO(3) existiert
eine eindeutige Eulerwinkelzerlegung im Sinne der folgenden sukzessiven Anwendung von Rotationen:
R = R(α, β, γ) = Rz (α)Ry (β)Rz (γ),
mit 0 ≤ α, γ < 2π und 0 ≤ β ≤ π (vgl. Definition 1.4).
1.2 Wigner-D-Funktionen
In diesem Abschnitt werden spezielle Funktionen betrachtet, welche von der Definitionsmenge SO(3) in die Zielmenge C der komplexen Zahlen abbilden. Diese Funktionen sind
also von der Form
f : SO(3) → C.
In Abschnitt 1.1, Satz 1.5, wurde mittels der Eulerwinkelzerlegung bereits ein Koordinatensystem definiert, in welchem mit solchen Funktionen gearbeitet werden kann. Über dem
Körper C ist das Tripel ({f : SO(3) → C}, +, 0) ein Vektorraum. Von besonderem Interesse sind im Folgenden die Wigner-D-Funktionen, welche eine Teilmenge der Funktionen
f : SO(3) → C bilden. Um diese zu definieren, folgt zunächst die Definition der Wigner-dFunktionen, mit deren Hilfe dann die Wigner-D-Funktionen definiert werden.
Definition 1.6 ([Kostelec und Rockmore, 2008, S. 156]). Die Wigner-d-Funktionen sind mit
(α,β)
Hilfe der Jacobi-Polynome Pn
m, m0
: [−1, 1] → R [Edmonds, 1996, S. 57ﬀ.] für l ∈ N0 und
= −l, . . . , l definiert als
s
dlm,m0 (β) := ξm,m0
(2l + 1)s!(s + u + v)!
2(s + u)!(s + v)!
β u
β v (u,v)
sin
cos
Ps (cos β),
2
2
mit 0 ≤ β ≤ π, wobei u := |m − m0 |, v := |m + m0 |, s := l − (u + v)/2 und
ξm,m0

1
:=
(−1)m0 −m
wenn m0 ≥ m,
wenn m0 < m
gilt.
Definition 1.7 ([Kostelec und Rockmore, 2008, S. 155]). Die Wigner-D-Funktionen sind
mittels der Wigner-d-Funktionen definiert als
l
Dm,m
0 (α, β, γ) :=
1 −imα l
0
e
dm,m0 (β)e−im γ ,
2π
mit l ∈ N0 und |m|, |m0 | ≤ l.
0 ≤ β ≤ π, 0 ≤ α, γ < 2π,
6
Mathematische Theorie
l
0
Für die Wigner-D-Funktionen gilt {Dm,m
0 | l ∈ N0 , |m|, |m | ≤ l} ⊂ {f : SO(3) → C},
d.h. diese Funktionen bilden von der Rotationsgruppe in die komplexen Zahlen ab und können daher in dem in Satz 1.5 für diese Funktionen definierten Koordinatensystem verwendet werden. Nun wird mittels der Wigner-D-Funktionen ein Untervektorraum von {f :
SO(3) → C} erzeugt.
Definition 1.8. Die lineare Hülle der Menge aller Wigner-D-Funktionen bis zu einem gegebenen Grad B spannt über dem Körper C einen Vektorraum HB auf, welcher definiert ist
als
n
o
l
0
HB := span Dm,m
0 l < B, |m|, |m | ≤ l .
Die Menge HB ist per Definition ein Untervektorraum von {f : SO(3) → C}. Es stellt sich
heraus, dass die Wigner-D-Funktionen bis zum Grad B eine Orthonormalbasis (ONB) von
HB bilden. Um dies zu begründen, wird der Begriﬀ der Länge und des Winkels für diese
Funktionen und somit ein Skalarprodukt der Form
h·, ·i : HB × HB → C
und die dazugehörige Norm k · k =
p
h·, ·i benötigt. Ein solches Skalarprodukt ist wie folgt
definiert.
Definition 1.9 ([Kostelec und Rockmore, 2008, S. 148]). Es seien f, g ∈ HB zwei beliebige
Funktionen. Das Skalarprodukt dieser beiden Funktionen ist definiert als
Z2π Zπ Z2π
(1)
hf, gi :=
f (α, β, γ)g(α, β, γ) sin β dγ dβ dα,
0
0
0
wobei g das komplex Konjugierte von g bezeichnet.
Der Vektorraum HB ist zusammen mit dem in Definition 1.9 definierten Skalarprodukt ein
Prähilbertraum. Weiterhin erfüllen die Wigner-D-Funktionen folgende Orthonormalitätsrelation.
l
Satz 1.10 ([Kostelec und Rockmore, 2008, S. 149]). Die Menge der Funktionen {Dm,m
0 |
l < B, |m|, |m0 | ≤ l} bildet ein System von orthonormalen Funktionen bezüglich des Skalarproduktes (1), d.h.
(2)
D
E
l1
l2
Dm
,
D
= δl1 ,l2 δm1 ,m2 δm01 ,m02
0
0
m2 ,m
1 ,m
1
2
für l1 , l2 < B ∈ N und |m1 |, |m01 | ≤ l1 sowie |m2 |, |m02 | ≤ l2 , wobei δ das Kronecker-DeltaSymbol bezeichnet.
Wigner-D-Funktionen
7
Somit bilden die Wigner-D-Funktionen bis zum Grad B eine ONB für HB . Diese Eigenschaft ist wichtig, da die Orthonormalität dieser Funktion es nun ermöglicht, die eindeutige
Fourier-Darstellung einer Funktion f ∈ HB zu bestimmen, indem die Skalarprodukte der
l
0
Funktion f mit den Elementen der Basis {Dm,m
0 | l < B, |m|, |m | ≤ l}, die sogenannten
Fourier-Koeﬃzienten von f , berechnet werden (siehe auch Abschnitt 1.3).
Abschließend werden zusätzliche Eigenschaften der Wigner-d-Funktionen genannt.
Die Jacobi-Polynome gehören zu den Orthogonalpolynomen. Eine Eigenschaft, welche diese
Orthogonalpolynome besitzen, ist, dass jedes Orthogonalpolynom eine Drei-Term-Rekursion erfüllt. Da die Wigner-d-Funktionen die Jacobi-Polynome als Faktor enthalten (vgl. Definition 1.6), erfüllen auch diese Funktionen eine Drei-Term-Rekursion.
Satz 1.11 ([Kostelec und Rockmore, 2008, S. 157]). Die Wigner-d-Funktionen erfüllen die
Drei-Term-Rekursion
dl+1
m,m0 (β)
mm0
=p
cos β −
dlm,m0 (β)
l(l + 1)
((l + 1)2 − m2 )((l + 1)2 − m02 )
(l + 1)(2l + 1)
p
(l + 1) (l2 − m2 )(l2 − m02 )
− p
dl−1
m,m0 (β)
2
2
2
02
l ((l + 1) − m )((l + 1) − m )
für l > 1, |m|, |m0 | ≤ l, wobei die folgenden Basisfälle die Rekursion initialisieren:
dll,m (β)
dl−l,m (β)
dlm,l (β)
dlm,−l (β)
s
(2l + 1)(2l)!
2(l + m)!(l − m)!
β l+m
β l−m
cos
− sin
,
2
2
s
(2l + 1)(2l)!
2(l + m)!(l − m)!
β l−m
β l+m
cos
sin
,
2
2
s
(2l + 1)(2l)!
2(l + m)!(l − m)!
β l+m
β l−m
cos
sin
,
2
2
s
(2l + 1)(2l)!
2(l + m)!(l − m)!
β l−m
β l+m
cos
− sin
.
2
2
=
=
=
=
Zusätzlich weisen diese Funktionen auch besondere Symmetrieeigenschaften auf [Edmonds,
1996, S. 59f.]. Diese Symmetrien sind gegeben durch
0
dlm,m0 (β) = (−1)m−m dl−m,−m0 (β)
m−m0
= (−1)
0
dlm0 ,m (β)
= (−1)l−m dl−m,m0 (π − β)
(3)
8
Mathematische Theorie
= (−1)l+m dlm,−m0 (π − β)
0
= (−1)l−m dl−m0 ,m (π − β)
= (−1)l+m dlm0 ,−m (π − β)
= dl−m0 ,−m (β)
für l ∈ N0 und |m|, |m0 | ≤ l.
1.3 Diskrete Fourier-Transformation auf SO(3) (DSOFT) und diskrete inverse Transformation (iDSOFT)
Es sei B ∈ N gegeben. Der in Abschnitt 1.2 definierte Vektorraum HB und das dazugehörige
Skalarprodukt h·, ·i ermöglichen es nun, für Funktionen h ∈ HB eine eindeutige FourierDarstellung anzugeben, welche durch
(4)
h=
B−1
X
l
l
X
X
l
ĥlm,m0 Dm,m
0
l=0 m=−l m0 =−l
gegeben ist [Koecher, 1997, S. 159]. Die Menge {ĥlm,m0 | l < B, |m|, |m0 | ≤ l} sind die
Fourier-Koeﬃzienten der Funktion h. Unter Nutzung der Eigenschaften des Skalarprodukts
und der Orthonormalität der Wigner-D-Funktionen folgt sofort
D
E
l
h, Dm,m
=
0
*B−1 j
j
X X X
+
j
l
ĥjk,k0 Dk,k
0 , Dm,m0
j=0 k=−j k0 =−j
(5)
=
B−1
X
j
j
X
X
D
E
j
l
ĥjk,k0 Dk,k
= ĥlm,m0
0 , Dm,m0
j=0 k=−j k0 =−j
für l < B, |m|, |m0 | ≤ l. In (5) wird die Orthonormalität der Wigner-D-Funktionen gel
0
nutzt, wodurch das Skalarprodukt einer Funktion aus {Dm,m
0 | l < B, |m|, |m | ≤ l} mit
sich selbst 1 und 0 sonst ergibt. Das heißt, es wird mittels des Skalarproduktes von h und den
Wigner-D-Funktionen, welche eine ONB des HB sind, eine Zerlegung der Funktion h bestimmt, sodass die Linearkombination der durch die entsprechenden Fourier-Koeﬃzienten
skalierten Wigner-D-Funktionen die Funktion h erzeugt. Aus der Orthornormalitätsrelation (2) folgt für alle l ≥ B, |m|, |m0 | ≤ l, sofort
ĥlm,m0 = 0.
Diskrete Fourier-Transformation auf SO(3) (DSOFT) und diskrete inverse
9
Transformation (iDSOFT)
Daher nennt man Funktionen h ∈ HB auch bandbegrenzt mit der Bandbegrenzung B (vgl.
[Kostelec und Rockmore, 2008, S. 151]). Aus (1) und (5) folgt nun:
Folgerung 1.12 (vgl. [Kostelec und Rockmore, 2008, S. 149]). Es sei h ∈ HB . Dann sind
die Fourier-Koeﬃzienten {ĥlm,m0 | l < B, |m|, |m0 | ≤ l} gegeben durch
ĥlm,m0
Z2π Zπ Z2π
=
0
0
l
h(α, β, γ)Dm,m
0 (α, β, γ) sin β dγ dβ dα
0
für l < B, |m|, |m0 | ≤ l.
Die bijektive Abbildung von h nach {ĥlm,m0 | l < B, |m|, |m0 | ≤ l} wird im Folgenden als
Fourier-Transformation auf SO(3) (SOFT) bzw. in Rückrichtung als inverse Fourier-Transformation auf SO(3) (iSOFT) bezeichnet.
Um die Integrale aus Folgerung 1.12 mittels eines Computers zu berechnen, müssen diese
durch endliche Summen ersetzt, d.h. diskretisiert werden. Daher muss die Funktion h auf
einem endlichen Gitter abgetastet werden. Im Bezug auf Folgerung 1.12 bedeutet dies, dass
die endlichen Summen, welche die Integrale ersetzen, die Form
X
X
X
l
0
0 0
w(α0 , β 0 , γ 0 )h(α0 , β 0 , γ 0 )Dm,m
0 (α , β , γ ),
α0 ∈Xα β 0 ∈Xβ γ 0 ∈Xγ
haben, wobei w eine Gewichtungsfunktion ist und X = Xα × Xβ × Xγ ⊂ SO(3) ein dreidimensionales Gitter bezeichnet, bei dem Xα , Xβ , Xγ die Abtastpunkte im Gitter entlang der
α- bzw. β-, bzw. γ-Achse sind (vgl. Satz 1.5).
Für Funktionen h ∈ HB kann nun eine Gewichtungsfunktion w mit einem Gitter X für
die exakte Berechnung der Fourier-Koeﬃzienten von h definiert werden. Entlang der α-, βund γ-Achse wird die Funktion h in äquidistanten Schritten abgetastet:
Definition 1.13 ([Kostelec und Rockmore, 2008, S. 151]). Es sei h ∈ HB . Dann ist die
Menge der Abtastpunkte des Gitters X = XB definiert als
XB :=
2πj1 π(2k + 1) 2πj2
R
,
,
=
2B
4B
2B
Rz
2πj1
π(2k + 1)
2πj2
Ry
Rz
2B
4B
2B
0 ≤ j1 , j2 , k < 2B ⊂ SO(3)
(vgl. Satz 1.4 und 1.5).
Die zu XB zugehörige Gewichtungsfunktion w = wB kann nun wie folgt definiert werden:
10
Mathematische Theorie
Definition 1.14 ([Kostelec und Rockmore, 2008, S. 151]). Die Gewichtungsfunktion wB ist
definiert als
2
w(R(j1 , k, j2 )) = wB (k) := sin
B
π(2k + 1)
4B
B−1
X
j=0
1
sin
2j + 1
π(2j + 1)(2k + 1)
4B
für R(j1 , k, j2 ) ∈ XB .
Diese Teile können nun zusammengesetzt werden, um die Fourier-Koeﬃzienten {ĥlm,m0 |
l < B, |m|, |m0 | ≤ l} aus Folgerung 1.12 mittels endlicher Summen exakt zu bestimmen.
Theorem 1.15 ([Kostelec und Rockmore, 2008, S. 151]). Es sei h ∈ HB . Dann können die
Fourier-Koeﬃzienten {ĥlm,m0 | l < B, |m|, |m0 | ≤ l} mittels
ĥlm,m0
(6)
=
π 2 2B−1
X 2B−1
X 2B−1
X
B
j1 =0 j2 =0 k=0
2πj1 π(2k + 1) 2πj2
wB (k) h
,
,
2B
4B
2B
×
l
Dm,m
0
2πj1 π(2k + 1) 2πj2
,
,
2B
4B
2B
für l < B, |m|, |m0 | ≤ l, berechnet werden.
Die Berechnung der Fourier-Koeﬃzienten mittels der endlichen Summen in Theorem 1.15
wird im Folgenden als diskrete Fourier-Transformation auf SO(3) (DSOFT) bezeichnet. Die
Rekonstruktion der Funktion h an den entsprechenden Stützstellen aus den Fourier-Koeﬃzienten {ĥlm,m0 | l < B, |m|, |m0 | ≤ l} mittels (4) wird dementsprechend als inverse
diskrete Fourier-Transformation auf SO(3) (iDSOFT) bezeichnet.
1.4 Schnelle Fourier-Transformation auf SO(3) (SOFFT) und schnelle inverse Transformation (iSOFFT)
Im Folgenden sei wieder B ∈ N gegeben. Eine mit B bandbegrenzte Funktion h ∈ HB
besitzt B(4B 2 −1)/3 Fourier-Koeﬃzienten [Kostelec und Rockmore, 2008, S. 164]. Die drei
Summen in (6) benötigen O(B 3 ) Operationen um einen dieser Koeﬃzienten zu berechnen.
Setzt man diese Teile zusammen, so benötigt man bei einer naiven Berechnung der FourierKoeﬃzienten mittels der DSOFT insgesamt O(B 6 ) Operationen.
Durch Umsortieren der Summanden kann die Anzahl der benötigten Operationen verringert werden. Diese Technik wird in [Kostelec und Rockmore, 2008, S. 152] als Separation der
Variablen bezeichnet. Zunächst werden dabei die Wigner-D-Funktionen durch ihre Darstellung von Exponential-Termen und der Wigner-d-Funktion ersetzt (vgl. Definition 1.7).
Schnelle Fourier-Transformation auf SO(3) (SOFFT) und schnelle inverse
11
Transformation (iSOFFT)
Anschließend wird die Reihenfolge der Summanden verändert. Gleichung (6) hat dann die
Form:
ĥlm,m0
2B−1
2B−1
2B−1
π X
π(2k + 1) X im0 2πj2 X im 2πj1
l
2B
=
w
(k)d
e
e 2B
0
B
m,m
(2B)2
4B
j2 =0
k=0
j1 =0
2πj1 π(2k + 1) 2πj2
×h
,
,
,
2B
4B
2B
mit l < B, |m|, |m0 | ≤ l. Eine vereinfachte Darstellung dieser Summen ist wie folgt gegeben:
2B−1 1 X
2πj1 π(2k + 1) 2πj2 im 2πj1
S1 (k, m, j2 ) :=
h
,
,
e 2B ,
2B
2B
4B
2B
|m| < B (7)
j1 =0
2B−1
1 X
0 2πj2
S2 (k, m, m ) :=
S1 (k, m, j2 )eim 2B ,
2B
0
|m|, |m0 | < B (8)
j2 =0
ĥlm,m0 = π
2B−1
X
wB (k)dlm,m0
k=0
π(2k + 1)
S2 (k, m, m0 ), |m|, |m0 | ≤ l < B. (9)
4B
Die Summen S1 und S2 können nun für festes k und j2 , bzw. festes k und m, mit jeweils
O(B log B) Operationen mittels einer schnellen Fourier-Transformation (FFT) berechnet
werden [Cormen et al., 2009, S. 901]. Durch zweifache Anwendung der inversen schnellen
Fourier-Transformation (iFFT) kann aus den Summen S2 die Funktion h an den entsprechenden Stützstellen mit O(B log B) Operationen rekonstruiert werden [Cormen et al.,
2009, S. 913]. Die Berechnung der Summen (9) wird wie folgt bezeichnet.
Definition 1.16 ([Kostelec und Rockmore, 2008, S. 153]). Für gegebene Ordnungen m und
m0 mit |m|, |m0 | < B ist die diskrete Wigner-d-Transformation (DWT) definiert als die Berechnung von
[v̂m,m0 ]l := π
2B−1
X
wB (k)dlm,m0
k=0
π(2k + 1)
[v]k
4B
(10)
0
für max{|m|, |m0 |} ≤ l < B, wobei v̂ ∈ C2B−max{|m|,|m |} der transformierte Vektor des
Eingangsvektors v ∈ C2B ist.
Für gegebene m und m0 benötigt die DWT O(B 2 ) Operationen zur Berechnung aller ĥlm,m0 ,
l = max{|m|, |m0 |}, . . . , B − 1. Zusammen mit den oben genannten FFTs ergibt dies
nun eine Anzahl von insgesamt O(B 4 ) Operationen für die Berechnung aller {ĥlm,m0 |
12
Mathematische Theorie
l < B, |m|, |m0 | ≤ l}. Im Folgenden wird diese Methode zur Berechnung der FourierKoeﬃzienten als schnelle Fourier-Transformation auf SO(3) (SOFFT) bezeichnet.
Sind die Ordnungen m und m0 gegeben, so lässt sich (9) auch als Matrix-Vektor-Produkt einer Matrix Wm,m0 ·G, welche die benötigten gewichteten Wigner-d-Funktionswerte enthält,
und einem Datenvektor v schreiben.
Satz 1.17 ([Kostelec und Rockmore, 2008, S. 153]). Es seien m und m0 mit |m|, |m0 | < B
gegeben, Wm,m0 eine Matrix, welche die Funktionswerte der Wigner-d-Funktionen enthält:
π(2k + 1)
0
Wm,m0 := djm,m0 cos
∈ R(B−max{|m|,|m |})×2B ,
j = max{|m|,|m0 |},...,B−1
4B
k = 0,...,2B−1
G eine Diagonalmatrix, auf deren Diagonale die Werte der Gewichtungsfunktion wB stehen:
G := diag(wB (0), . . . , wB (2B − 1)) ∈ R2B×2B
und v ∈ C2B ein gegebener Eingangsvektor. Dann ist der transformierte Vektor v̂m,m0 ∈
0
CB−max{|m|,|m |} der DWT gegeben durch das Matrix-Vektor-Produkt
v̂ = W · G · v.
Der rekonstruierte Eingangsvektor v eines gegebenen transformierten Vektors v̂m,m0 ist dann
gegeben durch die inverse diskrete Wigner-d-Transformation (iDWT) als das Matrix-VektorProdukt
v = W T · v̂.
Bei der Berechnung der SOFFT enthält der Vektor v̂ die Fourier-Koeﬃzienten {ĥlm,m0 | l =
max{|m|, |m0 |}, . . . , B − 1} der Funktion h,
(11)
max{|m|,|m0 |}+1 max{|m|,|m0 |} T
B−2
v̂ = ĥB−1
,
ĥ
,
.
.
.
,
ĥ
,
ĥ
,
0
0
0
0
m,m
m,m
m,m
m,m
für gegebene Ordnungen m, m0 mit |m|, |m0 | < B.
Mittels der iDWT und der iFFT (siehe oben) kann die Funktion h an den entsprechenden
Stützstellen aus den Fourier-Koeﬃzienten {ĥlm,m0 | l < B, |m|, |m0 | ≤ l} ebenfalls in
O(B 4 ) Operationen rekonstruiert werden. Dieses Verfahren wird im Folgenden als inverse
schnelle Fourier-Transformation auf SO(3) (iSOFFT) bezeichnet.
13
2
Implementierung und Parallelisierung
In diesem Kapitel wird die C++-Implementierung der SOFFT und iSOFFT vorgestellt. Dazu
werden in Abschnitt 2.1 die implementierten Algorithmen zunächst in Form von Pseudocode dargestellt und zur mathematischen Theorie in Zusammenhang gebracht. Abschnitt 2.2
gibt eine kurze Einführung in Open Multi-Processing (OpenMP), eine Spracherweiterung
der Programmiersprachen C, C++ und Fortran zur Parallelisierung von sequentiell implementierten Algorithmen. Anschließend wird in Abschnitt 2.3 die Forster-Methodik eingeführt. Die in Abschnitt 2.2 beschriebenen Konzepte werden dann zusammen mit der ForsterMethodik in Abschnitt 2.4 genutzt, um eine Parallelisierung der Implementierung aus Abschnitt 2.1 vorzunehmen. Abschließend wird in Abschnitt 2.5 die Zeitkomplexität und der
Arbeitsaufwand untersucht und die Optimalität der Parallelisierung nachgewiesen.
2.1 Sequentielle SOFFT und iSOFFT
Die sequentielle SOFFT sowie die sequentielle iSOFFT sind in Form von Pseudocode als Algorithmus 2.1 bzw. 2.2 dargestellt. Ein Algorithmus wird dann als sequentiell bezeichnet,
wenn die Anweisungen, aus welchen sich dieser Algorithmus zusammensetzt, in der vorgegebenen Reihenfolge nacheinander (in einer Sequenz) ausgeführt werden.
Es seien nun B ∈ N und h ∈ HB gegeben. Zunächst wird eine genaue Beschreibung der
Datentypen gegeben, welche bei der Implementierung der SOFFT/iSOFFT und der Darstellung als Pseudocode verwendet wurden. Mit Hilfe der Abtastpunkte aus Definition 1.13 wird
ein Gitter XB definiert, auf welchem die Funktion h abgetastet wird. Dieses Gitter bildet einen geometrischen Würfel, welcher in jeder der drei Dimensionen (α, β, γ) jeweils 2B Abtastpunkte enthält. Im folgenden Pseudocode enthält das dreidimensionale Datenarray G die
Funktionswerte von h über dem Gitter XB . Da bei der Implementierung die Ebenen senkrecht zur β-Achse und die Vektoren entlang der β-Achse in G benötigt werden, wird durch
die Notation Gk die k-te Ebene senkrecht zur β-Achse und durch Gα0 ,γ 0 der Vektor entlang
der β-Achse an der Stelle α0 in α-Richtung und γ 0 in γ-Richtung bezeichnet. Die Notation
G(n) bedeutet, dass G ein leeres dreidimensionales Array der Größe n×n×n ist. Neben dem
Datenarray G gibt es noch eine weitere Datenstruktur ĥ, in der die Fourier-Koeﬃzienten von
h gespeichert werden.
Nun wird der Algorithmus 2.1 genauer untersucht. Entsprechend (7), (8) und (9) erfolgt die
Berechnung der Fourier-Koeﬃzienten durch die SOFFT mehrstufig. Zunächst werden in den
Zeilen 10 bis 12 die Summen (7) und (8) für alle k = 0, . . . , 2B − 1 gemäß (9) berechnet.
Das heißt, alle Ebenen Gk senkrecht zur β-Achse werden nacheinander in α- und γ-Richtung
mittels der FFT transformiert. Die Indizes auf der α- und γ-Achse in G entsprechen nun den
14
Implementierung und Parallelisierung
1 /*
2
* Parameter:
3
*
in
G - Das Datenarray mit Funktionswerten einer Funktion
4
*
in
B - Bandbreite
5
*
out
ĥ - Das Datenarray mit den Fourier-Koeffizienten von
6
*/
h ∈ HB
h
7 function ĥ = SOFFT(G(2B), B) {
β -Achse
γ -Achse
8
// Transformiere jede Ebene Gk der
9
// und anschließend entlang der
10
α-Achse
for k = 0, . . . , 2B-1 seq do {
Transformiere Gk mittels FFT in
11
12
in G mittels FFT entlang der
α-Richtung
und anschließend in
γ -Richtung.
}
13
14
// Berechnung der Fourier-Koeffizienten für
15
Berechne aus G0,0 mittels der DWT
16
ĥl0,0
l = 0, . . . , B-1
17
für
18
in ĥ.
m = m0 = 0
und speichere diese mit den bereits berechneten Fourier-Koeffizienten
19
20
for m = 1, . . . , B-1 seq do {
// Berechne alle Fourier-Koeffizienten, bei welchen eine der Ordnungen gleich 0 ist,
21
22
// oder bei welchen der Betrag beider Ordnungen gleich ist.
23
Berechne aus
Gm,0 , G0,m , G−m,0 , G0,−m , Gm,m , G−m,−m , Gm,−m , G−m,m
24
mittels der DWT jeweils
25
26
ĥlm,0 , ĥl0,m , ĥl−m,0 , ĥl0,−m , ĥlm,m , ĥl−m,−m , ĥlm,−m , ĥl−m,m
l = m, . . . , B-1 und speichere diese mit den bereits
27
für
28
zienten in ĥ.
berechneten Fourier-Koeffi-
29
0
for m =
30
1, . . . , m
{
31
// Berechne übrige Fourier-Koeffizienten
32
Berechne aus
Gm,m0 , Gm0 ,m , G−m,−m0 , G−m0 ,−m , Gm0 ,−m , Gm,−m0 , G−m0 ,m , G−m,m0
33
mittels der DWT jeweils
34
35
ĥlm,m0 , ĥlm0 ,m , ĥl−m,−m0 , ĥl−m0 ,−m , ĥlm0 ,−m , ĥlm,−m0 , ĥl−m0 ,m , ĥl−m,m0
l = max{|m|, |m0 |}, . . . , B-1 und speichere diese mit den bereits
36
für
37
Fourier-Koeffizienten in ĥ.
}
38
39
berechneten
}
40 }
Algorithmus 2.1: Sequentielle SOFFT
Ordnungen m und m0 . Die Ordnung m stellt den Index in α-Richtung und die Ordnung m0
den Index in γ-Richtung dar. Anschließend wird aus G jeweils ein Vektor Gm,m0 entnommen
und mittels der DWT (10) wie in Satz 1.17 beschrieben in einen Vektor transformiert, welcher wie in (11) die Fourier-Koeﬃzienten von h enthält. Ist eine der beiden Ordnungen m
oder m0 negativ, so wird jeweils von 2B − 1 rückwärts gezählt.
Sequentielle SOFFT und iSOFFT
15
1 /*
2
* Parameter:
3
*
in
B - Bandbreite
4
*
in
ĥ - Das Datenarray mit den Fourier-Koeffizienten einer Funktion
5
*
out
G - Das Datenarray mit Funktionswerten der Funktion
6
*/
h ∈ HB
h
7 function G = ISOFFT(ĥ, B) {
8
// Ein leeres Datenarray für die Funktionswerte von
9
G = G(2B)
h
erstellen
10
11
// Berechnung des Datenvektors für
12
l
Berechne aus ĥ0,0 ,
l = 0, . . . , B-1,
m = m0 = 0
mittels iDWT den Datenvektor G0,0 .
13
14
for m = 1, . . . , B-1 seq do {
15
// Berechne alle Datenvektoren Gm,m0 , bei welchen eine der Ordnungen
16
// 0 ist, oder bei denen der Betrag beider Ordnungen gleich ist.
17
Berechne aus
l
l
l
l
l
l
l
m, m0
gleich
l
ĥm,0 , ĥ0,m , ĥ−m,0 , ĥ0,−m , ĥm,m , ĥ−m,−m , ĥm,−m , ĥ−m,m ,
18
l = m, . . . , B-1, mittels iDWT die entsprechenden Datenvektoren
Gm,0 , G0,m , G−m,0 , G0,−m , Gm,m , G−m,−m , Gm,−m , G−m,m .
19
20
21
0
for m = 1, . . . , m {
22
23
// Berechne übrige Datenvektoren Gm,m0
24
Berechne aus
l
l
l
l
l
l
0
l
l = max{|m|, |m |}, . . . , B-1, mittels iDWT die entsprechenden
Gm,m0 , Gm0 ,m , G−m,−m0 , G−m0 ,−m , Gm0 ,−m , Gm,−m0 , G−m0 ,m , G−m,m0 .
26
27
Datenvektoren
}
28
29
l
ĥm,m0 , ĥm0 ,m , ĥ−m,−m0 , ĥ−m0 ,−m , ĥm0 ,−m , ĥm,−m0 , ĥ−m0 ,m , ĥ−m,m0 ,
25
}
30
31
// Transformiere jeder Ebene Gk der
32
// anschließend entlang der
33
for k = 0, . . . , 2B-1 seq do {
Transformiere Gk mittels iFFT in
34
35
β -Achse
in G mittels iFFT entlang der
α-Achse
und
γ -Achse
α-Richtung
und anschließend in
γ -Richtung.
}
36 }
Algorithmus 2.2: Sequentielle iSOFFT
Zur Berechnung der Fourier-Koeﬃzienten von h werden die Symmetrien der Wigner-dFunktionen (3) genutzt, um die Anzahl der Auswertungen der Wigner-d-Funktionen in der
Implementierung zu reduzieren. Dazu laufen die For-Schleifen der Zeilen 20 und 30 jeweils
nur die positiven Wertebereiche von m und m0 ab. Die For-Schleife in Zeile 20 berechnet
zunächst alle Fourier-Koeﬃzienten, bei denen eine der beiden Ordnungen m und m0 gleich
0 ist, oder bei denen der Betrag beider Ordnungen gleich ist. Direkt im Anschluss werden
in der For-Schleife der Zeile 30 die übrigen Fourier-Koeﬃzienten berechnet, bis auf diejenigen mit den Ordnungen m = m0 = 0. Aus praktischen Gründen werden diese FourierKoeﬃzienten in Zeile 15 separat berechnet.
16
Implementierung und Parallelisierung
Die iSOFFT kann direkt aus der SOFFT abgeleitet werden und ist in Form von Pseudocode
als Algorithmus 2.2 dargestellt. Die DWT wird durch die in Satz 1.17 beschriebene iDWT
und die FFT durch die iFFT ersetzt. Zusätzlich wird die Reihenfolge aller Transformationen umgekehrt. Das heißt, es werden aus
ĥ
in den Zeilen 14 bis 29 zunächst die Fourier-
Koeﬃzienten von h entnommen und mittels der iDWT die entsprechenden FFT-transformierten Funktionswerte von h in den Ebenen
Gk
berechnet. Anschließend werden in den
Zeilen 33 bis 35 mittels der iFFT entlang der α- und anschließend entlang der γ-Achse in
den Ebenen Gk die ursprünglichen Funktionswerte von h in G rekonstruiert.
Um die FFT bzw. iFFT aus [Cormen et al., 2009, S. 901 bzw. S. 913] zu verwenden, wird stets
eine Zweierpotenz als Anzahl der Stützstellen benötigt. Dies impliziert insbesondere, dass
die Fourier-Koeﬃzienten einer Funktion h nur dann berechnet werden können, wenn die
Bandbreite B dieser Funktion ebenfalls eine Zweierpotenz ist. Um diese Einschränkung aufzuheben, wurde in der Implementierung auf die FFTW-Library [Frigo und Johnson, 2012]
zurückgegriﬀen. Die FFTW-Library stellt eine FFT-Version mit einer Zeitkomplexität von
O(n log n) für eine beliebige Anzahl n von Stützstellen zur Verfügung [Frigo und Johnson,
2012, S. 1].
2.2 Open Multi-Processing (OpenMP)
OpenMP ist eine Spracherweiterung der Programmiersprachen C, C++ und Fortran, mittels
welcher Codefragmente vom Compiler parallelisiert werden können [Quinn, 2003, S. 404].
Der Programmierer hat die Möglichkeit, mit Hilfe von Präprozessoranweisungen die entsprechend zu parallelisierenden Codefragmente zu markieren. Die Parallelisierung erfolgt im
Gegensatz zu anderen Programmierschnittstellen auf Ebene von Schleifen, während die Parallelisierung durch beispielsweise das Message Passing Interface (MPI) [Quinn, 2003, S. 95]
auf Prozessebene erfolgt, wodurch die Berechnungen ganzer Prozesse mittels Nachrichtenaustausch gesteuert werden müssen. Das heißt, jedes parallele Codefragment wird in einem
eigenen Prozess gekapselt und ausgeführt. Sofern ein Prozess Daten eines anderen Prozesses
benötigt, so wird ein Nachrichtenaustausch zwischen diesen Prozessen durchgeführt.
Das Konzept, welches OpenMP zugrunde liegt, wird als das Fork-Join-Modell bezeichnet [Cormen et al., 2009]. Dabei wird ein Programm über den gesamten Ausführungszeitraum betrachtet. Zu jedem Zeitpunkt der Ausführung kann sich die Anzahl der parallel auszuführenden Aktivitäten ändern. Abbildung 2.3 zeigt eine graphische Darstellung dieses Modells. Ein
Programm startet in einem Master-Thread und wird sequentiell abgearbeitet, bis ein Fork
eingeleitet wird. Eine vorher festgelegte Anzahl von Threads, das so genannte Team, wird erstellt, und jeder Thread erhält ein Arbeitspaket, welches er bearbeitet. Sobald ein Thread mit
seiner Arbeit fertig ist, erhält er entweder ein weiteres Arbeitspaket, oder verschmilzt mit
dem Master-Thread (Join), sofern keine weiteren Arbeitspakete berechnet werden müssen.
Entwurf paralleler Algorithmen mittels Forster-Methodik
Sequentieller Teil
Start
Erster paralleler Teil
Zweiter paralleler Teil
Innerhalb eines parallelen Pfades finden die
Operationen stets sequentiell statt. Es gibt
keinen Fork in einem Fork.
Ein neuer Fork kann nur
nach einem Join aller
Teammitglieder erfolgen.
Rechenoperation
(int a = b + 4;)
17
Ende
Abbildung 2.3: Beispiel der Parallelisierung eines Programmablaufs, dargestellt als Fork-Join-Modell.
Die Knoten stellen einzelne Berechnungen innerhalb des Programmablaufs dar. Der Pfad, welcher
durch die dicken Pfeile gebildet wird, stellt den Master-Thread dar. Parallele Codesegmente werden
dort eingeleitet, wo sich der Master-Thread in mehrere Ausführungspfade aufgliedert (Fork) und endet
dort, wo sich alle diese Pfade wieder mit dem Master-Thread vereinigt haben (Join). Jeder Pfad eines
Forks stellt die sequentielle Berechnung eines Threads dar.
Ein Arbeitspaket besteht aus einer Reihe von Berechnungen, welche sequentiell abgearbeitet
werden. Der Master-Thread wartet solange mit weiteren Berechnungen, bis alle Threads des
Teams ihre Arbeit erledigt haben. Eine Kommunikation zwischen den Threads ist nur durch
geteilte Speicherbereiche möglich, weshalb dieses Parallelisierungskonzept nur auf SharedMemory-Architekturen genutzt werden kann. In Shared-Memory-Architekturen verwenden
alle Prozessoren denselben Speicher.
2.3 Entwurf paralleler Algorithmen mittels Forster-Methodik
Bevor eine Parallelisierung der in Abschnitt 2.1 vorgestellten Algorithmen vorgenommen
wird, ist es nötig, einen parallelen Algorithmenentwurf anzufertigen. Zu diesem Zweck wird
die so genannte PCAM-Methodik, die in der Literatur häufig auch Forster-Methodik (nach
ihrem Erfinder Ian Forster) genannt wird, verwendet. Diese Entwurfsmethodik gliedert sich
in die vier Phasen Partitioning, Communication, Agglomeration und Mapping, welche entsprechend ihrer Reihenfolge durchlaufen werden (siehe Abbildung 2.4). Der so entstandene
parallele Algorithmus sollte in der Praxis mit einer geringeren Laufzeit berechnet werden
können als sein sequentielles Gegenstück. Aufgrund der Eigenschaften eines Mehrprozessorsystems kann dies allerdings nicht bei allen Algorithmen gewährleistet werden. Die genauen
Gründe werden später genannt (siehe auch Kapitel 4). Bevor eine detaillierte Beschreibung
der einzelnen Phasen erfolgt, wird zunächst der Begriﬀ des Overheads erklärt.
Definition 2.1 ([JáJá, 1992, Quinn, 2003, Roosta, 2000]). Die einzelnen, sequentiell arbeitenden Rechenwerke eines Parallelrechners werden im Folgenden als Einheiten bezeichnet.
Einheiten müssen bei der Berechnung paralleler Algorithmen explizit oder implizit miteinander kommunizieren. Diese Kommunikation kostet zusätzliche Zeit, die der entsprechende
18
Implementierung und Parallelisierung
sequentielle Algorithmus zur Berechnung des ihm zugrunde liegenden Problems nicht benötigt. Die zusätzlich benötigte Zeit, welche durch Kommunikation und zusätzliche Berechnungen anfällt, wird daher als Overhead bezeichnet. Das Ziel einer guten Parallelisierung
besteht darin, diesen Overhead so gering wie möglich zu halten. Zu diesem Zweck wurde
die Forster-Methodik entwickelt.
Die erste der vier Phasen wird als Partitioning bezeichnet. In dieser Phase wird versucht, so
viele Daten und Berechnungen wie möglich zu finden, die sich parallelisieren lassen. Diese Daten und Berechnungen werden dabei in kleine Arbeitspakete partitioniert. In OpenMP
kann die gewünschte Partitionierung durch die Präprozessoranweisungen #pragma
omp par-
allel for und #pragma omp parallel sections [The OpenMP Architecture Review Board,
2013, S. 96ﬀ.] angezeigt werden.
Im Anschluss an die Partitioning-Phase erfolgt die Communication-Phase. Nach der Identifizierung der Arbeitspakete für ein gegebenes Problem muss eine entsprechende Kommunikation zum Datenaustausch zwischen den Threads, die diese Arbeitspakete bearbeiten, bestimmt werden. Ein direkter Nachrichtenaustausch zwischen Threads ist im OpenMP-Standard nicht vorgesehen, sodass die einzige Möglichkeit der Kommunikation in der Nutzung
von geteilten Variablen (shared variables) besteht. Die Direktiven
private, lastprivate
shared, private, first-
und threadprivate [The OpenMP Architecture Review Board, 2013,
S. 44ﬀ.] ermöglichen es dem Compiler, Hinweise darauf zu geben, dass zwischen Threads
über die angegebenen Variablen kommuniziert wird.
Nachdem die Kommunikationswege zwischen den Arbeitspaketen bestimmt wurden, erfolgt
die Agglomeration (Ballung). Hier werden mehrere kleinere Arbeitspakete zu größeren Paketen zusammengefasst, um den Verwaltungsaufwand zur Bearbeitung der Arbeitspakete zu
verringern. Je besser die Agglomeration von Arbeitspaketen ist, desto geringer ist der Kommunikationsaufwand der größeren Pakete nach der Ballung. Die geringere explizite Kommunikation reduziert den Overhead, was wiederum zu einer geringeren Laufzeit des parallelisierten Algorithmus führt. Bei einer guten Agglomeration sind die resultierenden Pakete ähnlich groß. In OpenMP können dem Compiler Hinweise zur Agglomeration z.B. mittels der Direktiven if
(n > 100)
oder num_threads(8) [The OpenMP Architecture Review
Board, 2013, S. 48] gegeben werden. Eine weitere Möglichkeit, die Agglomeration einer Implementierung zu steuern, besteht darin, kleinere Schleifen sequentiell zu berechnen, oder
die Reihenfolge von Schleifen zu ändern, sodass agglomerierte Arbeitspakete direkt nacheinander abgearbeitet werden.
Die letzte Phase des parallelen Entwurfs ist das Mapping. Bei der Ausführung des Programms
werden die agglomerierten Arbeitspakete freien Einheiten zugeordnet. Optimalerweise sollten alle Einheiten zu jedem Zeitpunkt gleichmäßig ausgelastet sein. Die Umsetzung des Mappings erfolgt in OpenMP z.B. mittels der Anweisung
schedule
mit den Parametern
dyna-
Parallelisierung der SOFFT/iSOFFT mittels Forster-Methodik und OpenMP
19
Nebenläufigkeit und Skalierbarkeit
Communication
ng
tioni
Parti
Problem
Agglomeration
Mapping
Lokalität und Performanz
Abbildung 2.4: Schematische Darstellung der PCAM-Methodik nach [Quinn, 2003, S. 65]. Beginnend
mit einem Problem wird zunächst eine Partitionierung der Problemstellung in kleinere Arbeitspakete
vorgenommen. Im Anschluss wird eine Kommunikation zum Datenaustausch zwischen den Arbeitspaketen definiert. Im vorletzten Schritt werden die verbundenen Arbeitspakete zu größeren Arbeitspaketen geballt, sodass im letzten Schritt mittels des Mappings diese geballten Arbeitspakete auf die
Recheneinheiten verteilt werden können.
mic, static, guided, auto
und runtime [The OpenMP Architecture Review Board, 2013, S.
55ﬀ.].
Während in den Phasen Partitioning und Communication das Hauptaugenmerk darauf liegt,
möglichst viel Parallelität im gegebenen Problem zu erkennen, liegt der Fokus der letzten beiden Phasen, Agglomeration und Mapping, darauf, die theoretischen Überlegungen an die gegebenen Betriebsmittel anzupassen. Unter Betriebsmitteln versteht man in diesem Kontext
das System, auf welchem der Algorithmus ausgeführt werden soll. Je besser die Partitionierung ist und je weniger explizite Kommunikation notwendig ist, desto besser lässt sich ein
Problem parallelisieren und skalieren. In den Phasen der Agglomeration und des Mappings
kann durch geschickte Ballung der Arbeitspakete und der optimalen Ausnutzung der vorhandenen freien Einheiten dafür gesorgt werden, dass sich die Kommunikation zwischen
Paketen reduziert und sich die Laufzeit des parallelisierten Algorithmus verbessert. Gelingt
es beispielsweise, alle kleinen Arbeitspakete, die bei Berechnungen Abhängigkeiten zueinander haben (bei denen die bearbeitenden Threads miteinander kommunizieren müssen),
in ein großes Paket zu ballen, so muss nicht mehr explizit kommuniziert werden. Dies wird
auch als Lokalität bezeichnet.
2.4 Parallelisierung der SOFFT/iSOFFT mittels Forster-Methodik und
OpenMP
Alle Erläuterungen dieses Abschnitts beziehen sich auf die SOFFT aus Algorithmus 2.1 bzw.
iSOFFT aus Algorithmus 2.2. Gemäß Forster-Methodik wird in der Partitioning-Phase mit
20
Implementierung und Parallelisierung
m0
1
B−1
1
Thread 1 mit i = 1, . . . , (B − 2)/n − 1
Thread 2 mit i = (B − 2)/n, . . . , 2(B − 2)/n − 1
m
..
.
2
Thread n mit i = (B − 2) /n, . . . , (B − 1)(B − 2)/n − 1
B−1
Abbildung 2.5: Schematische Darstellung der Parallelisierung der äußeren For-Schleife in Algorithmus 2.1, Zeile 20. Der hellblau hinterlegte Bereich markiert die durch beide For-Schleifen abgedeckten
Indizes. Die roten Markierungen in m0 -Richtung gliedert die Wertebereiche des Index m in gleichmäßige Bereiche. Jeder Thread erhält bei Beginn der Berechnung einen dieser Bereiche und arbeitet alle
Indizes m und m0 innerhalb dieses Bereichs ab. Sofern dieser Thread mit seiner Arbeit fertig ist und
noch unbearbeitete Bereiche vorhanden sind, erhält er einen weiteren Bereich, welchen er abarbeitet.
Die rot markierten Indexpaarungen auf der Diagonalen sind bereits berechnete Indizes.
m0
1
B−1
1
0
T1 : Thread 1 mit j = 1, . . . , m /n − 1
0
0
T2 : Thread 2 mit j = m /n, . . . , 2m /n − 1
.
.
.
m
0
0
Tn : Thread n mit j = (B − 2)m /n, . . . , (B − 1)m /n − 1
B−1
…
T1
T2
Tn
Abbildung 2.6: Schematische Darstellung der Parallelisierung der inneren For-Schleife in Zeile 30 aus
Algorithmus 2.1. Der hellblau hinterlegte Bereich markiert die durch beide For-Schleifen abgedeckten
Indizes. Die roten Markierungen in m-Richtung gliedert die Wertebereiche des Index m0 in gleichmäßige Bereiche. Jeder Thread erhält bei Beginn der Berechnung einen dieser Bereiche und arbeitet alle
Indizes m und m0 innerhalb dieses Bereichs ab. Sofern dieser Thread mit seiner Arbeit fertig ist und
noch unbearbeitete Bereiche vorhanden sind, erhält er einen weiteren Bereich, welchen er abarbeiten
kann. Die rot markierten Indexpaarungen auf der Diagonalen sind bereits berechnete Indizes.
dem Entwurf einer Parallelisierung begonnen. Da OpenMP, wie bereits in Abschnitt 2.2 erklärt, auf Ebene von For-Schleifen parallelisiert, müssen eben diese im Algorithmus zunächst
identifiziert werden. Neben den For-Schleifen in den Zeilen 20 und 30 besteht die Möglichkeit, die DWT, welche durch ein Matrix-Vektor-Produkt beschrieben wird, nebenläufig zu
berechnen. Bei der Parallelisierung wurde die DWT aus praktischen Gründen nicht berücksichtigt. Die Gründe dafür werden in Kapitel 4 erläutert. Somit fällt die Entscheidung darauf,
die Berechnungen der For-Schleifen in den Zeilen 20 und 30 parallel auszuführen. Um dies
Parallelisierung der SOFFT/iSOFFT mittels Forster-Methodik und OpenMP
21
zu bewerkstelligen, werden nun zwei naheliegende Parallelisierungsansätze vorgestellt. Beide
dieser Parallelisierungsstrategien weisen Schwachstellen auf, die im späteren Verlauf durch
eine Kombination beider Ansätze beseitigt werden. Da in OpenMP keine Hinweise für eine
Parallelisierung in einem parallel ausgeführten Codesegment gegeben werden können ( Es
”
gibt keinen Fork in einem Fork“, siehe Abbildung 2.3), muss somit eine Entscheidung getroﬀen werden, welche der beiden For-Schleifen parallelisiert werden soll.
Die erste Strategie besteht darin, die äußere For-Schleife (Algorithmus 2.1, Zeile 20) zu parallelisieren. Dies bedeutet, es wird ein Team von Threads initialisiert, auf welches der Indexwertebereich von m gleichmäßig aufgeteilt wird. Die Anweisungen innerhalb der ForSchleife werden für die entsprechenden Indizes sequentiell berechnet. Schematisch entspricht
dies Abbildung 2.5. Es ist gut zu sehen, dass der entsprechend hellblau markierte Workload
pro Thread unregelmäßig verteilt ist. Die resultierende Parallelisierung kann daher die entsprechende Partitionierung nicht eﬃzient ausnutzen.
Die zweite Strategie verfolgt den Ansatz, die innere For-Schleife (Algorithmus 2.1, Zeile 30)
zu parallelisieren. Schematisch entspricht dies Abbildung 2.6. Zu sehen ist, dass nun jeder
Thread pro Iteration der äußeren For-Schleife gleichmäßig mit Arbeit versorgt wird, wodurch das Problem des vorherigen Ansatzes gelöst wäre. Diese Strategie hat allerdings einen
Nachteil. Durch die Parallelisierung der inneren For-Schleife werden die Berechnungen der
äußeren For-Schleife bei der Parallelisierung nicht berücksichtigt. Auch hier ist die Partitionierung gröber als möglich.
Es wird also eine feinere Partitionierung benötigt, welche dem Compiler mehr Spielraum
beim Scheduling der Arbeitspakete einräumt. Da ohne eine explizite Änderung am Algorithmus diese Möglichkeit nicht besteht, wird dazu zunächst eine Umformung des Algorithmus 2.1 vorgenommen. Dieser enthält nach der Umformung keine verschachtelten ForSchleifen mehr, wodurch die Arbeitspakete der Partitionierung wesentlich kleiner sind. Um
dies zu erreichen, wird im Folgenden ein geometrischer Lösungsansatz verfolgt, mittels welchem die äußere For-Schleife umgeformt wird.
Die dreieckige Form des Indexwertebereichs in den Abbildungen 2.5 und 2.6 indiziert in Algorithmus 2.1 eine direkte Abhängigkeit des Index m0 der For-Schleife in Zeile 30 zum Index
m der For-Schleife in Zeile 20. Um diese Abhängigkeit aufzulösen, wird die For-Schleife in
drei Schritten umgeformt. Zunächst wird die For-Schleife aus Algorithmus 2.1, Zeile 20 in
zwei For-Schleifen aufgeteilt. Die entsprechende Änderung ist in Algorithmus 2.7, Zeile 20
bzw. Zeile 31 zu sehen. Die For-Schleife in Zeile 20 enthält nun die Berechnungen für die
Fourier-Koeﬃzienten, welche vorher in den Zeilen 23 bis 28 in Algorithmus 2.1 berechnet
wurden, und die For-Schleife in Zeile 31 die innere For-Schleife aus Algorithmus 2.1, Zeile 30. Im zweiten Schritt wird die Abhängigkeit zwischen den For-Schleifen der Zeilen 31 und
32 in Algorithmus 2.7 aufgelöst. Zu diesem Zweck wird die dreieckige Form des Indexwerte-
22
Implementierung und Parallelisierung
19
for m = 1, . . . , B-1 seq do {
20
21
// Berechne alle Fourier-Koeffizienten, bei welchen eine der Ordnungen gleich 0 ist,
22
// oder bei welchen der Betrag beider Ordnungen gleich ist.
23
Berechne aus
Gm,0 , G0,m , G−m,0 , G0,−m , Gm,m , G−m,−m , Gm,−m , G−m,m
24
mittels der DWT jeweils
25
26
ĥlm,0 , ĥl0,m , ĥl−m,0 , ĥl0,−m , ĥlm,m , ĥl−m,−m , ĥlm,−m , ĥl−m,m
l = m, . . . , B-1 und speichere diese mit den bereits
27
für
28
zienten in ĥ.
berechneten Fourier-Koeffi-
}
29
30
for m = 1, . . . , B-1 seq do {
31
0
for m =
32
1, . . . , m
{
33
// Berechne übrige Fourier-Koeffizienten
34
Berechne aus
Gm,m0 , Gm0 ,m , G−m,−m0 , G−m0 ,−m , Gm0 ,−m , Gm,−m0 , G−m0 ,m , G−m,m0
35
mittels der DWT jeweils
36
37
ĥlm,m0 , ĥlm0 ,m , ĥl−m,−m0 , ĥl−m0 ,−m , ĥlm0 ,−m , ĥlm,−m0 , ĥl−m0 ,m , ĥl−m,m0
l = max{|m|, |m0 |}, . . . , B-1 und speichere diese mit den bereits
38
für
39
Fourier-Koeffizienten in ĥ.
berechneten
}
40
}
41
42 }
Algorithmus 2.7: Die letzten 24 Zeilen der ersten Umformung von Algorithmus 2.1. Die Zeilen 1 bis 19
sind identisch zu denen in Algorithmus 2.1. Die äußere For-Schleife aus Algorithmus 2.1, Zeile 20 wurde
in zwei For-Schleifen aufgeteilt, sodass die neue For-Schleife in Zeile 20 die Berechnungen der FourierKoeﬃzienten enthält und die zweite For-Schleife in Zeile 31 den Index m für die innere For-Schleife in
Zeile 32 zur Verfügung stellt.
bereichs in eine rechteckige Form überführt. Dies bedeutet, dass die For-Schleifen nach der
Transformation unabhängig voneinander sind. Die geometrische Veranschaulichung dieser
Überlegung aus Abbildung 2.8 spiegelt das Vorgehen wider. Der Anschaulichkeit halber wird
hierbei die Bandbreite B als ungerade angenommen. Andernfalls sind leichte Anpassungen
nötig.
Wird der Wertebereich von m in den Abbildungen 2.5 und 2.6 halbiert, so ist leicht zu sehen,
dass der untere Teil des gefärbten dreieckigen Wertebereichs der For-Schleifen-Indizes genau der ungefärbten Fläche in der oberen Hälfte des Quadrates entspricht. Daraus lässt sich
ableiten, dass die Anzahl der Indexpaarungen von m und m0 im hellblau gefärbten Bereich
genau so groß ist wie die Anzahl der Paarungen aller Werte der oberen Hälfte des gesamten Vierecks. Daher ist es möglich, zwei neue voneinander unabhängige Indizes i und j auf
die voneinander abhängigen Indizes m und m0 abzubilden. Für die Indizes i und j gilt nun
i = 1, . . . , d(B − 1)/2e, j = 1, . . . , B − 1. Mit Hilfe entsprechender Fallunterscheidungen
können nun aus allen i und j alle gültigen Paarungen von Werten für m und m0 rekonstruiert
werden. Die Umformung ist in Algorithmus 2.9 dargestellt.
Parallelisierung der SOFFT/iSOFFT mittels Forster-Methodik und OpenMP
23
3 j = 1, . . . , B − 1
m
1
3
B−1
i = 1, . . . ,
2
B−1
2
1
1
B−1
m0
2
B−1
1
B−1
2
−1
2
Abbildung 2.8: Geometrischer Ansatz der For-Schleifen-Optimierung. Die Punkte stellen die jeweiligen Kombinationen der Indizes m und m0 dar. 1 Der dreieckige Indexwertebereich wird zunächst in
m-Richtung halbiert. 2 Anschließend wird der Wertebereich der unteren Quadrathälfte um beide Achsen gedreht. 3 Der gespiegelte Bereich füllt nun den leeren Bereich in der oberen Quadrathälfte aus.
Der rechteckige Bereich, der die obere Quadrathälfte widerspiegelt, kann mittels zweier unabhängiger
Indizes i und j mit entsprechenden For-Schleifen abgelaufen werden. Die rot markierten Indexpaarungen auf der Diagonalen sind bereits berechnete Indizes und werden durch eine Indexverschiebung in
m0 -Richtung übersprungen.
In einem letzten Umformungsschritt werden nun die verschachtelten For-Schleifen mit den
Indizes i und j durch eine einzelne For-Schleife ersetzt. Dazu wird die exakte Anzahl der
Indexpaarungen für m und m0 benötigt. Aus der dreieckigen Form in den Abbildungen 2.5
und 2.6 lässt sich ableiten, dass die Anzahl der Indizes die Summe der ersten B − 2 Zahlen
bildet. Somit kann ein Index k = 0, . . . , (B − 1)(B − 2)/2 − 1 auf die benötigten Indizes i
und j und damit auch auf m und m0 abgebildet werden. Die Werte für i und j können nun
mittels einer einfachen Division bzw. des dazugehörigen Teilerrests ermittelt werden. Algorithmus 2.10 auf Seite 26 stellt den vollständig umgeformten Algorithmus dar. Da nun ausschließlich einfache For-Schleifen vorhanden sind (Algorithmus 2.10 Zeilen 10, 21 und 33),
können diese parallel berechnet werden. Somit ist die Partitionierung abgeschlossen. Die
Größe der entstandenen Arbeitspakete entspricht einer Indexpaarung bzw. im Fall der ForSchleife aus Zeile 10 einer Vektorlänge von Gm,m0 und ist daher nahezu minimal.
Im Anschluss an die Partitioning-Phase folgt nun die Communication-Phase. Die DWT benötigt außer den Ordnungen m und m0 , welche als jeweilige Indexpaarung in den Arbeitspaketen enthalten sind, keine weiteren Daten der anderen Berechnungen. Es ist folglich keine
24
Implementierung und Parallelisierung
30
for i =
31
1, . . . , d(B-1)/2e
seq do {
for j = 1, . . . , B-1 seq do {
32
33
// rekonstruiere
34
if (j > i) {
35
m = B - i
36
m = B - j
m
und
m0
0
} else {
37
38
m = i + 1
39
m = j
0
40
}
41
// Berechne übrige Fourier-Koeffizienten
42
Berechne aus
Gm,m0 , Gm0 ,m , G−m,−m0 , G−m0 ,−m , Gm0 ,−m , Gm,−m0 , G−m0 ,m , G−m,m0
43
mittels der DWT jeweils
44
45
ĥlm,m0 , ĥlm0 ,m , ĥl−m,−m0 , ĥl−m0 ,−m , ĥlm0 ,−m , ĥlm,−m0 , ĥl−m0 ,m , ĥl−m,m0
l = max{|m|, |m0 |}, . . . , B-1 und speichere diese mit den bereits
46
für
47
Fourier-Koeffizienten in ĥ.
berechneten
}
48
}
49
50 }
Algorithmus 2.9: Die letzten 21 Zeilen der zweiten Umformung von Algorithmus 2.1. Die Zeilen 1 bis
30 sind identisch zu denen in Algorithmus 2.7. Die beiden For-Schleifen in den Zeilen 31 und 32 aus
Algorithmus 2.7 wurden durch zwei unabhängige For-Schleifen ersetzt. Im Schleifenrumpf werden vor
der Berechnung der Fourier-Koeﬃzienten die Indizes m und m0 rekonstruiert.
explizite Kommunikation nötig. Die Datenstrukturen ĥ und G sind geteilte Datenstrukturen.
In der späteren Implementierung erfolgt der Zugriﬀ durch die entsprechenden Einheiten
exklusiv, d.h. jede Einheit liest und schreibt ausschließlich in eigenen Bereichen innerhalb
dieser Datenstrukturen. Es gibt also keine Überschneidung dieser Bereiche zwischen den
Einheiten und somit ist auch hier keine Kommunikation nötig. Algorithmen dieser Art werden als exclusive read exclusive write (EREW) bezeichnet [JáJá, 1992, S. 15].
Die vorletzte Phase des Entwurfs ist die Agglomeration-Phase. Da OpenMP für die Parallelisierung verwendet wird, besteht die einzige Möglichkeit der Agglomeration darin, die Parallelisierung erst ab einer bestimmten Problemgröße vorzunehmen. Weitere Erläuterungen
dazu folgen in Kapitel 4.
Zuletzt erfolgt die Mapping-Phase. Da die Arbeitspakete sehr klein sind, wird mit Hilfe von
OpenMP jede Iteration der For-Schleifen aus Algorithmus 2.10 in den Zeilen 10, 21 und 33
nebenläufig berechnet. Zusätzlich dazu wird dem Compiler ein dynamisches Scheduling angezeigt. Dazu wird die OpenMP-Direktive schedule(dynamic) [The OpenMP Architecture
Review Board, 2013, S. 57] genutzt.
Der parallele Algorithmus skaliert mit variierender Problemgröße und Anzahl der vorhandenen Einheiten sehr gut. Der Grund dafür ist die nahezu minimale Größe der Arbeitspakete und die nicht notwendige explizite Kommunikation. Die geringe Kommunikation führt
Analyse der Parallelisierung
25
wiederum zu einer hohen Lokalität der Daten. Die feine Granularität der Partitionierung in
Arbeitspakete macht es dem Compiler darüber hinaus sehr einfach, ein gutes Scheduling für
die Bearbeitung dieser Arbeitspakete zu finden.
Wie in Abschnitt 2.1 beschrieben, kann die iSOFFT direkt aus der SOFFT abgeleitet werden,
indem die DWTs durch iDWTs und FFTs durch iFFTs ersetzt werden und die Reihenfolge
der For-Schleifen umgekehrt wird. Dies trifft ebenfalls auf die parallelisierte SOFFT zu. Somit
folgt aus Algorithmus 2.10 direkt die parallelisierte iSOFFT, dargestellt in Algorithmus 2.11
auf Seite 27. Die Eigenschaften dieses Algorithmus bezüglich Nebenläufigkeit, Skalierbarkeit,
Lokalität und Performanz sind identisch zu denen der parallelisierten SOFFT.
2.5 Analyse der Parallelisierung
Um die Parallelisierung aus Abschnitt 2.4 zu analysieren, werden im Folgenden die Begriﬀe
des Speedup, der Eﬃzienz, der Laufzeit und der Optimalität paralleler Algorithmen eingeführt. Im Anschluss daran wird die Güte der Parallelisierung mit Hilfe einiger dieser Maße
bewertet. Der Speedup und die Eﬃzienz werden in Kapitel 4 genutzt.
Definition 2.2 ([JáJá, 1992, Quinn, 2003, Roosta, 2000]). Die Rechenzeit eines parallelen
Algorithmus A bei einer Eingabe x und genau p Einheiten wird mit TpA (x) bezeichnet. Die
Rechenzeit TpA (n) bezeichnet die maximale Rechenzeit von A über allen Eingaben x der
Länge n.
Da die genaue Anzahl der später verwendeten Einheiten häufig nicht bekannt ist, werden
Algorithmen in der Theorie unter der Annahme, eine unbegrenzte Anzahl von Einheiten sei
verfügbar, analysiert. Stehen für die Berechnung eines parallelen Algorithmus unbegrenzt
viele Einheiten zur Verfügung, so kann die Zeitkomplexität dieses Algorithmus nicht schlechter sein, als die Berechnung mit einer vorher festgelegten Anzahl von Einheiten.
Definition 2.3 ([JáJá, 1992, Quinn, 2003, Roosta, 2000]). Die Rechenzeit eines parallelen
Algorithmus A, der beliebig viele Einheiten verwenden darf, wird mit
A
T∞
(n) := inf TpA (n)
p∈N
bezeichnet.
Neben der parallelen Rechenzeit eines Algorithmus für ein gegebenes Problem wird im Folgenden die sequentielle Rechenzeit dieses Algorithmus als Maßstab genutzt.
Definition 2.4 ([JáJá, 1992]). Die Zeitkomplexität des schnellsten (in der Regel unbekannten) sequentiellen Algorithmus A zur Lösung eines Problems mit Eingaben der Länge n wird
im Folgenden mit T A,∗ (n) bezeichnet.
26
Implementierung und Parallelisierung
1 /*
2
* Parameter:
3
*
in
G - Das Datenarray mit Funktionswerten einer Funktion
4
*
in
B - Bandbreite
5
*
out
ĥ - Das Datenarray mit den Fourier-Koeffizienten von
6
*/
h ∈ HB
h
7 function ĥ = SOFFT(G(2B), B) {
β -Achse
γ -Achse
8
// Transformiere jede Ebene Gk der
9
// und anschließend entlang der
10
α-Achse
for k = 0, . . . , 2B-1 par do {
Transformiere Gk mittels FFT in
11
12
in G mittels FFT entlang der
α-Richtung
und anschließend in
γ -Richtung.
}
13
14
// Berechnung der Fourier-Koeffizienten für
15
Berechne aus G0,0 mittels der DWT
16
ĥl0,0
l = 0, . . . , B-1
17
für
18
in ĥ.
m = m0 = 0
und speichere diese mit den bereits berechneten Fourier-Koeffizienten
19
20
// Berechnung der Fourier-Koeffizienten mit $m=0$, $m'=0$ oder
21
for m = 1, . . . , B-1 par do {
|m| = |m0 |
// Berechne alle Fourier-Koeffizienten, bei welchen eine der Ordnungen gleich 0 ist,
22
23
// oder bei welchen der Betrag beider Ordnungen gleich ist.
24
Berechne aus
Gm,0 , G0,m , G−m,0 , G0,−m , Gm,m , G−m,−m , Gm,−m , G−m,m
25
mittels der DWT jeweils
26
27
ĥlm,0 , ĥl0,m , ĥl−m,0 , ĥl0,−m , ĥlm,m , ĥl−m,−m , ĥlm,−m , ĥl−m,m
l = m, . . . , B-1 und speichere diese mit den bereits
28
für
29
zienten in ĥ.
30
berechneten Fourier-Koeffi-
}
31
32
// Berechnung aller übrigen Fourier-Koeffizienten
33
for k = 0, . . . , (B-1)(B-2)/2-1 par do {
m
und
m0
34
// Rekonstruiere
35
if (k % (B-1) + 1 > k / (B-1) + 1) {
36
m = B - k / (B-1) + 1
37
m = B - k % (B-1) + 1
0
} else {
38
39
m = k / (B-1) + 2
40
m = k % (B-1) + 1
0
41
}
42
// Berechne übrige Fourier-Koeffizienten
43
Berechne aus
Gm,m0 , Gm0 ,m , G−m,−m0 , G−m0 ,−m , Gm0 ,−m , Gm,−m0 , G−m0 ,m , G−m,m0
44
mittels der DWT jeweils
45
46
ĥlm,m0 , ĥlm0 ,m , ĥl−m,−m0 , ĥl−m0 ,−m , ĥlm0 ,−m , ĥlm,−m0 , ĥl−m0 ,m , ĥl−m,m0
l = max{|m|, |m0 |}, . . . , B-1 und speichere diese mit den bereits
47
für
48
Fourier-Koeffizienten in ĥ.
49
}
50 }
Algorithmus 2.10: Parallelisierte SOFFT
berechneten
Analyse der Parallelisierung
1 /*
2
* Parameter:
3
*
in
4
*
in
ĥ - Das Datenarray mit den Fourier-Koeffizienten einer Funktion
5
*
out
G - Das Datenarray mit Funktionswerten der Funktion
6
*/
B - Bandbreite
h ∈ HB
h
7 function G = ISOFFT(ĥ, B) {
8
// Ein leeres Datenarray für die Funktionswerte von
9
G = G(2B)
h
erstellen
10
11
// Berechnung des Datenvektors für
12
l
Berechne aus ĥ0,0 ,
l = 0, . . . , B-1,
m = m0 = 0
mittels iDWT den Datenvektor G0,0 .
13
14
// Berechnung der Fourier-Koeffizienten mit $m=0$, $m'=0$ oder $|m|=|m'|$
15
for m = 1, . . . , B-1 par do {
16
// Berechne alle Datenvektoren Gm,m0 , bei welchen eine der Ordnungen
17
// 0 ist, oder bei denen der Betrag beider Ordnungen gleich ist.
18
Berechne aus
l
l
l
l
l
l
gleich
l
ĥm,0 , ĥ0,m , ĥ−m,0 , ĥ0,−m , ĥm,m , ĥ−m,−m , ĥm,−m , ĥ−m,m ,
19
l = m, . . . , B-1, mittels iDWT die entsprechenden Datenvektoren
Gm,0 , G0,m , G−m,0 , G0,−m , Gm,m , G−m,−m , Gm,−m , G−m,m .
20
21
22
l
m, m0
}
23
24
// Berechnung aller übrigen Fourier-Koeffizienten
25
for k = 0, . . . , (B-1)(B-2)/2-1 par do {
m
und
m0
26
// Rekonstruiere
27
if (k % (B-1) + 1 > k / (B-1) + 1) {
28
m = B - k / (B-1) + 1
29
m = B - k % (B-1) + 1
0
} else {
30
31
m = k / (B-1) + 2
32
m = k % (B-1) + 1
0
33
}
34
// Berechne übrige Datenvektoren Gm,m0
35
Berechne aus
l
l
l
l
l
l
l
0
l = max{|m|, |m |}, . . . , B-1, mittels iDWT die entsprechenden
Gm,m0 , Gm0 ,m , G−m,−m0 , G−m0 ,−m , Gm0 ,−m , Gm,−m0 , G−m0 ,m , G−m,m0 .
37
38
39
l
ĥm,m0 , ĥm0 ,m , ĥ−m,−m0 , ĥ−m0 ,−m , ĥm0 ,−m , ĥm,−m0 , ĥ−m0 ,m , ĥ−m,m0 ,
36
Datenvektoren
}
40
41
// Transformiere jeder Ebene Gk der
42
// anschließend entlang der
43
for k = 0, . . . , 2B-1 par do {
Transformiere Gk mittels iFFT in
44
45
β -Achse
in G mittels iFFT entlang der
α-Achse
γ -Achse
α-Richtung
und anschließend in
}
46 }
Algorithmus 2.11: Parallelisierte iSOFFT
γ -Richtung.
und
27
28
Implementierung und Parallelisierung
Wie bereits in [Kostelec und Rockmore, 2008, S. 153] beschrieben, besteht die Möglichkeit,
die DWT mittels einer schnellen diskreten Polynom-Transformation [Driscoll et al., 1997]
mit einer Zeitkomplexität von O(B log2 B) zu berechnen. Ersetzt man in der SOFFT die
DWT durch diese schnelle diskrete Polynom-Transformation, so sinkt die Zeitkomplexität des gesamten Algorithmus von O(B 4 ) auf O(B 3 log2 B). Somit hat diese modifizierte SOFFT die kleinste bekannte Zeitkomplexität zur Berechnung der Fourier-Koeﬃzienten
von Funktionen f ∈ HB . Da sich die Parallelisierung allerdings auf den in Abschnitt 1.4 beschriebenen SOFFT bzw. iSOFFT Algorithmus bezieht, wird im Folgenden T SOFFT,∗ (8B 3 ) ∈
O(B 4 ) für die sequentielle SOFFT/iSOFFT angenommen. Beachte, dass in diesem Kontext
die Eingabewörter eine Länge von 2B · 2B · 2B = 8B 3 haben.
Definition 2.5 ([JáJá, 1992]). Die Arbeit W A (x) eines parallelen Algorithmus A bei einer
Eingabe x ergibt sich aus der Summe der Anzahl von (unbegrenzt verfügbaren) Einheiten,
welche in jedem Zeitschritt tatsächlich arbeiten. Analog zu der Zeitkomplexität TpA (n) ist
die Arbeit W A (n) die maximale Arbeit über allen Eingaben x der Länge n.
Um die Güte eines parallelen Algorithmus zu bewerten, werden neben der Zeitkomplexität
und der Arbeit dieses Algorithmus die Maße des Speedup und der Eﬃzienz benötigt. Die
folgenden Definitionen dieser Maße basieren auf den Definitionen in [JáJá, 1992, Quinn,
2003, Roosta, 2000].
Definition 2.6. Der Speedup SpA (n) eines parallelen Algorithmus A ist definiert als
SpA (n) :=
T A,∗ (n)
,
TpA (n)
wobei n die Länge der Eingaben bezeichnet.
Definition 2.7. Die Eﬃzienz EpA (n) eines parallelen Algorithmus A ist definiert als
EpA (n) :=
SpA (n)
,
p
wobei n die Länge der Eingaben bezeichnet.
Damit später gezeigt werden kann, dass die Parallelisierung der SOFFT und somit auch der
iSOFFT aus Abschnitt 2.4 optimal ist, wird nun der Begriﬀ der Optimalität für parallele Algorithmen definiert.
Definition 2.8 ([JáJá, 1992, S. 32]). Ein paralleler Algorithmus A heißt optimal, falls für ihn
W A ∈ Θ T A,∗
Analyse der Parallelisierung
29
gilt, wobei Θ(f ) := O(f ) ∩ Ω(f ) mit
O(f ) := { g | ∃k > 0 ∃n0 ∈ N ∀n ≥ n0 : g(n) ≤ k · f (n) },
Ω(f ) := { g | ∃k > 0 ∃n0 ∈ N ∀n ≥ n0 : g(n) ≥ k · f (n) }
bezeichnet (f und g sind hierbei zwei Laufzeitfunktionen).
Anders als bei sequentiellen Algorithmen spielt die Zeitkomplexität TpA in Definition 2.8
keine Rolle beim Begriﬀ der Optimalität. Im Sequentiellen gilt ein Algorithmus erst dann
als optimal, wenn die Zeitkomplexität eine untere Schranke für das gegebene Problem ist. In
[JáJá, 1992] werden zwei Arten der Optimalität definiert. Um Missverständnisse durch Nutzung des Optimalitätsbegriﬀs zu vermeiden, wird die Optimalität aus Definition 2.8 auch als
Arbeitsoptimalität bezeichnet. Die Arbeit W SOFFT eines parallelen Algorithmus A für ein gegebenes Problem entspricht der Anzahl der Operationen, die alle Einheiten zusammen bei
Berechnung von A benötigen. Die Anzahl der Berechnungsschritte zur Lösung desselben
Problems mit Hilfe des schnellsten sequentiellen Algorithmus entspricht der Zeitkomplexität T A,∗ . Folglich wird unter Berücksichtigung von Definition 2.8 der parallele Algorithmus
A als arbeitsoptimal bezeichnet, wenn die Anzahl der benötigten Berechnungsschritte des
parallelen Algorithmus A identisch zu denen der schnellsten bekannten sequentiellen Implementierung ist.
Nun kann die parallelisierte SOFFT bzw. iSOFFT mit Hilfe dieser Maße analysiert werden.
Zunächst wird T SOFFT (8B 3 ) ∈ O(B 4 ) überprüft, wobei T SOFFT die Laufzeitfunktion der sequentiellen SOFFT bezeichnet. Die Korrektheit des Algorithmus muss nicht explizit gezeigt
werden, da diese sich aus der Korrektheit der mathematischen Definition ergibt.
In Algorithmus 2.1, Zeilen 10 bis 12 werden 2B ·4B FFTs berechnet. Der geometrische Würfel, welcher durch G dargestellt wird, enthält 2B zur β-Achse senkrechte Ebenen. Jede diese
Ebenen enthält jeweils 2B Werte in α- und γ-Richtung. Die entsprechende Ebene wird wie
in Abschnitt 2.1 beschrieben jeweils entlang der α- und γ-Achse mittels der FFT transformiert. Die Zeitkomplexität der FFT entlang der genannten Achsen beträgt O(2B log(2B))
Somit ergibt sich für diesen Abschnitt eine Zeitkomplexität von
2 · 2B · 2B · O(2B log(2B)) = O(16B 3 log(2B)) = O(B 3 log B).
In Algorithmus 2.1, Zeilen 23 bis 28 werden (B − 2) · 8 DWTs berechnet. Anschließend
werden die berechneten Fourier-Koeﬃzienten in ĥ gespeichert. Eine DWT ist das Produkt
einer Matrix und eines Vektors. Da die Wigner-d-Matrix Wm,m0 , welche für die Transformation genutzt wird, die Größe (B − max{|m|, |m0 |}) × 2B hat, ist diese für die Parameter m = m0 = 0 mit B × 2B am größten. Daher wird für die Laufzeitanalyse aller
30
Implementierung und Parallelisierung
DWTs die Matrix mit einer Größe von B × 2B abgeschätzt, wodurch eine Zeitkomplexität
von O(2B 2 ) = O(B 2 ) für eine DWT entsteht. Das Speichern der Fourier-Koeﬃzienten
kann mit einer linearen Zeitkomplexität in O(B) durchgeführt werden, da der FourierKoeﬃzienten-Vektor gemäß (11) für m = m0 = 0 genau B Werte enthält. Somit ergibt sich
für die Zeilen 23 bis 28 eine Zeitkomplexität von
(B − 1) · (8 · O(B 2 ) + 8 · O(B)) = (B − 1) · O(B 2 ) = O(B 3 ).
In Algorithmus 2.1, Zeilen 32 bis 37 werden 8·(B −2)(B −1)/2 DWTs berechnet. Anschließend werden die berechneten Fourier-Koeﬃzienten in ĥ gespeichert. Somit ist die Zeitkomplexität dieser Berechnungen
1
1
(B − 2)(B − 1) · (8 · O(B 2 ) + 8 · O(B)) = (B − 2)(B − 1) · O(B 2 ) = O(B 4 ).
2
2
Um die Zeitkomplexität des gesamten Algorithmus 2.1 zu bestimmen, wird nun die Summe
der Komplexitäten aller For-Schleifen zu der Zeitkomplexität der einzelnen DWT in Zeile 15
addiert. Dies ergibt nun eine Zeitkomplexität von
T SOFFT,∗ (8B 3 ) ∈ O(B 4 ) + O(B 3 log B) + O(B 3 ) + O(B 2 ) = O(B 4 ).
Die berechnete Zeitkomplexität stimmt damit mit der in [Kostelec und Rockmore, 2008]
beschriebenen Komplexität überein und wird im Folgenden benutzt, um die Arbeitsoptimalität von Algorithmus 2.10 nachzuweisen. Für die Zeitkomplexität T iSOFFT der iSOFFT gilt
nun T SOFFT = T iSOFFT , da wie bereits in Abschnitt 2.1 beschrieben die FFTs durch iFFTs und
die DWTs durch iDWTs ersetzt werden und im Anschluss die Reihenfolge der For-Schleifen
umgekehrt wird. Bei einer iDWT wird die Wigner-d-Matrix transponiert, daher ist die Zeitkomplexität identisch zu der einer DWT. Die Zeitkomplexität der iFFT entspricht der der
FFT.
Im nächsten Schritt wird die Zeitkomplexität von Algorithmus 2.10 bestimmt. Die Gesamtlaufzeit setzt sich wieder aus der Zeitkomplexität der drei For-Schleifen, sowie der Berechnung der Koeﬃzienten für m = m0 = 0 in Zeile 15 zusammen.
In Algorithmus 2.10, Zeilen 10 bis 12 werden alle Iteration der For-Schleife von 2B Einheiten gleichzeitig berechnet. Da zunächst eine der Achsen in α-Richtung oder γ-Richtung
mittels der FFT transformiert wird und im Anschluss die andere Achse, werden insgesamt
nur zwei Zeitschritte pro Iteration für die komplette Transformation der Ebene Gk benötigt.
Analyse der Parallelisierung
31
Da jede Ebene Gk , d.h. jede Iteration der For-Schleife, simultan berechnet wird, ergibt dies
eine Zeitkomplexität von
2 · O(2B log(2B)) = O(B log B).
In Algorithmus 2.10, Zeilen 21 bis 30 werden alle B − 2 Iterationen parallel berechnet. Die
Zeitkomplexität um ein Matrix-Vektor-Produkt parallel zu berechnen beträgt O(log n) [JáJá,
1992, S. 12]. In der Implementierung wurde auf eine Parallelisierung des Matrix-VektorProdukts verzichtet. Die Gründe dafür werden in Kapitel 3 genannt. Das Speichern der Fourier-Koeﬃzienten in ĥ kann mit B Einheiten in einem Zeitschritt erfolgen. Somit ergibt sich
für diese For-Schleife eine Laufzeit von
8 · O(log B) + 8 = O(log B).
In Algorithmus 2.10, Zeilen 33 bis 49 werden alle (B − 2)(B − 1)/2 Iterationen der ForSchleife parallel berechnet. Analog zur vorherigen For-Schleife werden dort acht DWTs berechnet. Im Anschluss daran werden die B Fourier-Koeﬃzienten in ĥ von B Einheiten in
einem Zeitschritt gespeichert. Dies resultiert ebenfalls in einer Zeitkomplexität von
8 · O(log B) + 8 = O(log B).
Insgesamt ergibt dies also zusammen mit der DWT in Algorithmus 2.10 Zeile 15 eine Zeitkomplexität von
SOFFT
T∞
(8B 3 ) ∈ O(B log B) + 2 · O(log B) + O(log B) = O(B log B).
SOFFT
Somit gilt T∞
(8B 3 ) ∈ O(B log B).
Die ermittelte Zeitkomplexität entspricht der einer sequentiellen FFT. Sofern ein entsprechender paralleler Algorithmus für die FFT verwendet wird, kann die Zeitkomplexität für die
parallele SOFFT auf O(log B) bei gleichbleibender Arbeit reduziert werden [JáJá, 1992, S.
384, Theorem 8.5]. In der Implementierung wurde allerdings auf Funktionen aus der FFTWLibrary zurückgegriﬀen. Diese Library stellt nebenläufige Funktionen zur Berechnung der
FFT zur Verfügung. Die Dokumentation liefert allerdings keine Hinweise zur parallelen Zeitkomplexität dieser Algorithmen, sodass für die obere Schranke die Zeitkomplexität einer
sequentiellen FFT genutzt wurde. Zuletzt wird gezeigt, dass die parallele SOFFT aus Algorithmus 2.10 arbeitsoptimal ist.
Theorem 2.9. Die parallelisierte SOFFT sowie die parallelisierte iSOFFT aus Algorithmus 2.10
bzw. Algorithmus 2.11 sind mit einer Arbeit von W SOFFT (8B 3 ) ∈ O(B 4 ) arbeitsoptimal.
32
Implementierung und Parallelisierung
Beweis. Analog zur Laufzeit der parallelen SOFFT wird die Arbeit W SOFFT bestimmt, indem
die Arbeit der drei For-Schleifen aus Algorithmus 2.10, Zeilen 10, 21 und 33 sowie die Arbeit
für die Berechnung der Fourier-Koeﬃzienten m = m0 = 0 in Zeile 15 aufsummiert wird.
Die For-Schleife in Zeile 10 durchläuft B Iterationen, in denen jeweils 2B Einheiten parallel
in zwei Zeitschritten die FFT entlang der α- und γ-Achsen berechnen. Dies ergibt eine Arbeit
von
B · 2B = 4B 2 = O(B 2 ).
In der For-Schleife aus Zeile 21 werden für die Berechnung eines Matrix-Vektor-Produkts
in logarithmischer Zeit B 2 Einheiten benötigt. Dies ergibt somit eine Arbeit von O(B 2 ) für
eine DWT. Die For-Schleife benötigt B − 1 Iterationen, welche von B − 1 Einheiten parallel berechnet werden. Für das Speichern der resultierenden Fourier-Koeﬃzienten werden
nochmals B Einheiten benötigt. Insgesamt ergibt sich für diese For-Schleife also eine Arbeit
von
(B − 1)(8 · O(B 2 ) + O(B)) = (B − 1) · O(B 2 ) = O(B 3 ).
Für die letze For-Schleife in Zeile 33 werden nun (B −1)(B −2)/2 Iterationen von verschiedenen Einheiten parallel berechnet. Die Arbeit der DWTs und das Speichern der FourierKoeﬃzienten beträgt weiterhin O(B 2 ) bzw. O(B). Dies führt zu einer Arbeit von
1
1
(B − 1)(B − 2)(8 · O(B 2 ) + O(B)) = (B 2 − 3B + 2) · O(B 2 ) = O(B 4 ).
2
2
Zuletzt muss noch die Summe über die Arbeit der einzelnen Teile sowie die der DWT für
die Ordnungen m = m0 = 0 in Zeile 15 gebildet werden:
W SOFFT (B) ∈ O(B 4 ) + O(B 3 ) + 2 · O(B 2 ) = O(B 4 ).
Es gilt somit W SOFFT (8B 3 ) ∈ O(B 4 ). Da W SOFFT ∈ Θ(T SOFFT,∗ ) mit T SOFFT,∗ ∈ O(B 4 ) gilt, ist
dieser Algorithmus optimal.
Wie bereits bei der Analyse der Zeitkomplexität beschrieben, ist die Zeit bzw. die Arbeit der
SOFFT identisch zu der Zeit bzw. der Arbeit der iSOFFT. Die Arbeitsoptimalität der SOFFT
impliziert somit die Arbeitsoptimalität der iSOFFT.
33
3
Numerische Beispiele
Die in Kapitel 2 vorgestellte Parallelisierung der SOFFT und iSOFFT wird nun mit Hilfe von
verschiedenen Benchmarks getestet. Zu diesem Zweck werden im Folgenden drei Benchmarks präsentiert. Zunächst wird ein Benchmark zur Analyse des Speedup der parallelisierten SOFFT bzw. iSOFFT vorgestellt. Anschließend wird ein Benchmark zur Analyse der
Eﬃzienz und zuletzt ein Benchmark zur Rechengenauigkeit der DWT aus Kapitel 1, Abschnitt 1.4, Satz 1.17 beschrieben. Abschließend werden in Abschnitt 3.2 die Ergebnisse dieser Benchmarks dargestellt.
3.1 Benchmark zu Speedup, Eﬃzienz und Rechengenauigkeit
Der Speedup-Benchmark misst den Speedup der parallelisierten SOFFT und iSOFFT (vgl.
Definition 2.6). Dazu misst der Benchmark die entsprechenden Rechenzeiten und berechnet daraus den Speedup. Benchmark 3.1 zeigt den Messvorgang. Bei der Ausführung dieses
Benchmarks werden die Bandbreiten 16, 32, 64 und 128 durchlaufen. Zuerst werden dabei
alle benötigten Datenstrukturen erzeugt. Das heißt, für den Messvorgang wird ein Gitter G
sowie die Datenstruktur ĥ angelegt. Für jede der Bandbreiten B wird zunächst eine Menge
an zufälligen Fourier-Koeﬃzienten gleichverteilt auf [−1, 1] × [−i, i] generiert. Der Pseudozufallszahlengenerator, welcher für diesen und alle weiteren Benchmarks verwendet wird,
ist der Mersenne-Twister MT19937 [Matsumoto und Nishimura, 1998]. Aus der Menge der
generierten Fourier-Koeﬃzienten werden zunächst mit Hilfe der iSOFFT im Gitter
G
die
Funktionswerte der entsprechenden Funktion h ∈ HB rekonstruiert. Im Anschluss daran
werden fünf sequentielle SOFFT durchgeführt, wobei die mittlere Laufzeit T SOFFT,∗ bestimmt
wird und als Referenz für die Berechnung des Speedup gespeichert wird. Dann wird mit Hilfe eines For-Schleifen-Index p die Anzahl der zu verwendenden Einheiten für die parallele
SOFFT gesteuert. Für jede Kombination aus Bandbreite B und Anzahl p an Einheiten wird
nun jeweils fünf Mal der SOFFT-Algorithmus ausgeführt. Danach wird der Mittelwert der
Laufzeit dieser fünf Durchgänge gebildet. Mit Hilfe der durchschnittlichen Laufzeit der parallelen SOFFT und T SOFFT,∗ kann nun der Speedup SpA (B) berechnet werden.
Um die Messungen für die iSOFFT durchzuführen, wird dieser Benchmark umformuliert
(siehe Benchmark 3.2). Genau wie beim Speedup-Benchmark für die SOFFT werden beim
Benchmark für die iSOFFT zunächst die benötigten Datenstrukturen instanziiert. Daraufhin
wird für die aktuelle Bandbreite B ∈ {16, 32, 64, 128} eine Menge an zufälligen FourierKoeﬃzienten generiert. Aus diesen Fourier-Koeﬃzienten werden mittels einer sequentiellen iSOFFT die Funktionswerte der entsprechenden Funktion h ∈ HB in
G
rekonstru-
iert. Die mittlere Laufzeit T iSOFFT,∗ für die sequentielle iSOFT wird aus fünf Durchgängen
34
Numerische Beispiele
(1) Erzeuge gleichverteilte Fourier-Koeﬃzienten ĥ aus [−1, 1] × [−i, i] für h ∈ HB .
(2) Rekonstruiere h in G aus ĥ mittels iSOFFT.
(3) Führe jeweils 5 sequentielle SOFFTs durch und bestimme die durchschnittliche Berechnungsdauer T SOFFT,∗ .
(4) Führe jeweils 5 parallele SOFFTs durch und bestimme die durchschnittliche Berechnungsdauer TpSOFFT .
(5) Wiederhole Schritt (4) für alle möglichen Anzahlen p an Einheiten, beginnend bei 2 und
endend bei der maximal möglichen Anzahl.
(6) Berechne für jede Anzahl p an Einheiten den dazugehörigen Speedup SpSOFFT (B) (vgl. Definition 2.6).
(7) Wiederhole Schritt (1) bis (6) für alle Bandbreiten B ∈ {16, 32, 64, 128}.
Benchmark 3.1: Benchmark zur Messung des Speedups. Die Messungen dieses Benchmarks beziehen
sich auf die parallelisierte SOFFT aus Algorithmus 2.10.
(1) Erzeuge gleichverteilte Fourier-Koeﬃzienten ĥ aus [−1, 1] × [−i, i] für h ∈ HB .
(2) Führe jeweils 5 sequentielle iSOFFTs durch und bestimme die durchschnittliche Berechnungsdauer T iSOFFT,∗ .
(3) Führe jeweils 5 parallele iSOFFTs durch und bestimme die durchschnittliche Berechnungsdauer TpiSOFFT .
(4) Wiederhole Schritt (3) für alle möglichen Anzahlen p an Einheiten, beginnend bei 2 und
endend bei der maximal möglichen Anzahl.
(5) Berechne für jede Anzahl p an Einheiten den dazugehörigen Speedup SpiSOFFT (B).
(6) Wiederhole Schritt (1) bis (5) für alle Bandbreiten B ∈ {16, 32, 64, 128}.
Benchmark 3.2: Benchmark zur Messung des Speedups. Die Messungen dieses Benchmarks beziehen
sich auf die parallelisierte iSOFFT aus Algorithmus 2.11.
der sequentiellen iSOFFT bestimmt. Anschließend wird mittels einer Schleife und einem
dazugehörigen Index p die Anzahl der Einheiten beginnend bei 2 bis zur maximal möglichen Anzahl durchlaufen. Analog zum Speedup-Benchmark für die SOFFT wird für jede
Paarung der Anzahl p an Einheiten und der aktuellen Bandbreite B fünf Mal der parallele iSOFFT-Algorithmus ausgeführt. Nach den fünf Durchläufen wird die durchschnittliche Berechnungsdauer der parallelen iSOFFT ermittelt. Mit Hilfe der durchschnittlichen
Berechnungsdauer der parallelisierten iSOFFT und T iSOFFT,∗ wird der Speedup SpA (B) berechnet.
Als nächstes wird der Eﬃzienz-Benchmark zur Messung der Eﬃzienz der Parallelisierung
von SOFFT und iSOFFT, dargestellt in den Benchmarks 3.3 und 3.4, vorgestellt. Es werden
wie beim Speedup-Benchmark Messungen für die Bandbreiten 16, 32, 64 und 128 durchgeführt. Für jede dieser Bandbreiten werden zuerst die nötigen Datenstrukturen instanziiert.
Anschließend werden gleichverteilte Fourier-Koeﬃzienten aus der Menge [−1, 1]×[−i, i] er-
Benchmark zu Speedup, Effizienz und Rechengenauigkeit
35
(1) Erzeuge gleichverteilte Fourier-Koeﬃzienten ĥ aus [−1, 1] × [−i, i] für h ∈ HB .
(2) Rekonstruiere h in G aus ĥ mittels iSOFFT.
(3) Führe fünf sequentielle SOFFTs durch und ermittle die Zeit T SOFFT,∗ .
(4) Führe jeweils 5 parallele SOFFTs durch und bestimme die durchschnittliche Berechnungsdauer TpSOFFT .
(5) Wiederhole Schritt (4) für alle möglichen Anzahlen p an Einheiten, beginnend bei 2 und
endend bei der maximal möglichen Anzahl.
(6) Berechne für jede Anzahl p an Einheiten die dazugehörige Eﬃzienz EpSOFFT (B) (vgl. Definition 2.7).
(7) Wiederhole Schritt (1) bis (6) für alle Bandbreiten B ∈ {16, 32, 64, 128}.
Benchmark 3.3: Benchmark zur Messung der Eﬃzienz. Die Messungen dieses Benchmarks beziehen
sich auf die parallelisierte SOFFT aus Algorithmus 2.10.
(1) Erzeuge gleichverteilte Fourier-Koeﬃzienten ĥ aus [−1, 1] × [−i, i] für h ∈ HB .
(2) Führe fünf sequentielle iSOFFTs durch und ermittle die Zeit T iSOFFT,∗ .
(3) Führe jeweils 5 parallele iSOFFTs durch und berechne die durchschnittliche Berechnungsdauer TpiSOFFT .
(4) Wiederhole Schritt (3) für alle möglichen Anzahlen p an Einheiten, beginnend bei 2 und
endend bei der maximal möglichen.
(5) Berechne für jede Anzahl p an Einheiten die dazugehörige Eﬃzienz EpiSOFFT (B).
(6) Wiederhole Schritt (1) bis (5) für alle Bandbreiten B ∈ {16, 32, 64, 128}.
Benchmark 3.4: Benchmark zur Messung der Eﬃzienz. Die Messungen dieses Benchmarks beziehen
sich auf die parallelisierte iSOFFT aus Algorithmus 2.11.
zeugt. Aus diesen Fourier-Koeﬃzienten werden mittels des iSOFFT-Algorithmus die Funktionswerte von h über dem Gitter
G
rekonstruiert. Mittels der rekonstruierten Funktions-
werte wird im Anschluss daran fünf mal die sequentielle SOFFT ausgeführt und die mittlere Ausführungszeit T SOFFT,∗ bestimmt. Danach wird mittels einer For-Schleife und dem
dazugehörigen Schleifen-Index p die Anzahl der möglichen Einheiten für die parallelisierte
SOFFT durchlaufen. Innerhalb dieser Schleife werden nun fünf Messungen der Rechenzeit
TpSOFFT (B), welche die parallelisierte SOFFT für jedes mögliche p und die aktuelle Bandbreite
B ∈ {16, 32, 64, 128} benötigt, durchgeführt. Dabei wird auch hier eine Variable genutzt,
welche die Summe der Ausführungszeiten enthält. Nachdem die fünf Messungen vorgenommen wurden, wird die durchschnittliche Ausführungszeit bestimmt. Abschließend wird die
Eﬃzienz EpA (B) gemäß Definition 2.7 berechnet. Der Benchmark der iSOFFT ist bis auf wenige Änderungen identisch zu dem der SOFFT. Zuerst werden die benötigten Datenstrukturen erzeugt. In Anschluss daran werden gleichverteilte Fourier-Koeﬃzienten aus der Menge
[−1, 1] × [−i, i] generiert. Diese Fourier-Koeﬃzienten werden nun wieder genutzt, um mit-
36
Numerische Beispiele
(1) Erzeuge Gewichtsmatrix G und Wigner-d-Matrix Wm,m0 für gegebenes B, m und m0 .
(2) Erzeuge gleichverteilte Fourier-Koeﬃzienten v̂ für h ∈ HB (vgl. (11)).
T
(3) Rekonstruiere Datenvektor s aus v̂ mittels iDWT durch s = Wm,m
0 · v̂.
(4) Berechne ĝ mittels DWT durch ĝ = Wm,m0 · G · s.
(5) Berechne Diﬀerenzvektor d = ĝ − v̂.
(6) Berechne den absoluten Fehler ∆d und relativen Fehler δd .
(7) Führe die Schritte (2) bis (6) 1.000 mal durch und berechne den mittleren absoluten Fehler
∆d und den mittleren relativen Fehler δd .
(8) Führe Schritt (7) für die Bandbreiten 16, 32, 64, 128 und die Ordnungen m = m0 = 0,
m = B/2 mit m0 = 0, m = m0 = B/2 durch.
Benchmark 3.5: Benchmark zur Messung der Rechengenauigkeit der DWT aus Satz 1.17.
tels der sequentiellen iSOFFT die Funktionswerte der entsprechenden Funktion h über dem
Gitter G zu rekonstruieren. Die mittlere Zeit T iSOFFT,∗ wird durch fünf sequentiellen SOFFTs
bestimmt. Danach wird wieder mittels einer Schleife und einem Index p jede mögliche Anzahl an Einheiten durchlaufen. In jeder Iteration werden fünf Messungen der Rechenzeit des
parallelisierten iSOFFT-Algorithmus für die aktuelle Bandbreite B durchgeführt. Es wird anschließend der Mittelwert der Rechenzeit über diese Durchläufe gebildet. Daraus wird dann
die Eﬃzienz EpA (B) berechnet.
Abschließend wird der Genauigkeits-Benchmark zur Messung der Genauigkeit der DWT aus
Satz 1.17 vorgestellt (siehe Benchmark 3.5). Dieser Benchmark misst den Rechenfehler, welcher auftritt, nachdem ein Datenvektor durch die DWT in Fourier-Koeﬃzienten und anschließend durch eine iDWT wieder in den ursprünglichen Datenvektor transformiert wurde. Die Messungen erfolgen für vorgegebene Ordnungen m und m0 und Bandbreiten 16,
32, 64 und 128. Zunächst wird für die aktuelle Bandbreite B die Gewichtsmatrix G aufgestellt. Anschließend wird die Wigner-d-Matrix Wm,m0 erstellt. Diese Wigner-d-Matrix wird
mit der Gewichtsmatrix gewichtet. Zudem wird eine Kopie der ungewichteten Wigner-dMatrix für die iDWT transponiert. Anschließend wird für jede Bandbreite jeweils 1.000
mal ein Vektor v̂ mit zufälligen Fourier-Koeﬃzienten (vgl. (11)) zuerst in einen Datenvektor und anschließend zurück in einen Fourier-Koeﬃzienten-Vektor transformiert. Jedes mal
wird die Diﬀerenz des ursprünglichen Fourier-Koeﬃzienten-Vektors und des rekonstruierten Fourier-Koeﬃzienten-Vektors berechnet. Aus diesem Diﬀerenzvektor d wird der absolute Fehler ∆d durch die euklidische Norm
∆d := kdk2 =
q
[d]21 + [d]22 + · · · + [d]2B−max{|m|,|m0 |}−1
des Diﬀerenzvektors berechnet. Neben dem absoluten Fehler ∆d wird auch der relative Fehler δd berechnet. Der relative Fehler berechnet sich aus der euklidischen Norm des Diﬀeren-
Ergebnisse der Benchmarks
37
zenvektors d geteilt durch die euklidische Norm des ursprünglichen Fourier-KoeﬃzientenVektors,
q
[d]21 + [d]22 + · · · + [d]2B−max{|m|,|m0 |}−1
∆d
δd :=
=q
.
2 + · · · + [v̂]2
kv̂k2
[v̂]2
+
[v̂
]
0
1
2
B−max{|m|,|m |}−1
Diese Fehler werden über alle 1.000 Iterationen gemittelt und durch ∆d als δd gespeichert.
3.2 Ergebnisse der Benchmarks
In diesem Abschnitt sind die Ergebnisse der im vorherigen Abschnitt vorgestellten Benchmarks dargestellt. Diese Ergebnisse wurden mit einer Intel Core i7-4470 Quad-Core CPU
mit einer Taktung von 3,4 GHz erzeugt (siehe auch [The Intel Corporation, 2015]). Mittels
Hyperthreading kann diese CPU Berechnungen mit acht Einheiten simultan durchführen.
Der benutzte Rechner verfügt über 8 GB Hauptspeicher. Das System wurde mit der Linux
Distribution Ubuntu mit der Version 14.03.3 LTS (64 Bit) betrieben. Zum Kompilieren wurde der g++-Compiler mit der Versionsnummer 4.8 und OpenMP-Unterstützung verwendet,
sowie der C++-Standard 11 mit der Optimierungsstufe 3 [The GCC-Team, 2015]. In allen
Benchmarks wurde die sequentielle DWT und iDWT verwendet. Darüber hinaus wurde die
SOFFT und iSOFFT erst ab einer Bandbreite von B = 20 parallel ausgeführt. Die Gründe
dafür werden in Kapitel 4 näher erläutert.
Abbildung 3.6 zeigt die Ergebnisse des Speedup-Benchmarks. Entlang der x-Achse sind die
8
8
7
7
6
6
Speedup
Speedup
Anzahl der Einheiten aufgetragen. Entlang y-Achse sind die zugehörigen Speedup-Werte
5
4
5
4
3
3
2
2
1
1
1
2
3
4
5
Einheiten
6
7
8
1
2
3
4
5
6
7
8
Einheiten
Abbildung 3.6: Ergebnisse des Speedup-Benchmarks. Im linken Diagramm sind die Ergebnisse der
SOFFT und im rechten Diagramm die Ergebnisse der iSOFFT zu sehen. Die orangene Kurve ( ) stellt
den Speedup für die Bandbreite 16 dar. Die weiteren Kurven stellen den Speedup für die Bandbreite
32 ( ), die Bandbreite 64 ( ) und die Bandbreite 128 ( ) dar. Die schwarz gepunktet Linie zeigt den
optimalen Speedup.
Numerische Beispiele
1
1
0.8
0.8
0.6
0.6
Eﬃzienz
Eﬃzienz
38
0.4
0.4
0.2
0.2
0
0
1
2
3
4
5
6
7
8
1
2
3
Einheiten
4
5
6
7
8
Einheiten
Abbildung 3.7: Ergebnisse des Eﬃzienz-Benchmarks. Im linken Diagramm ist die Eﬃzienz der SOFFT
zu sehen. Analog zeigt das rechte Diagramm die Eﬃzienz der iSOFFT. Die orangene Kurve ( ) stellt
die Eﬃzienz für die Bandbreite 16 dar. Die weiteren Kurven stellen die Eﬃzienz für die Bandbreite 32
( ), die Bandbreite 64 ( ) und die Bandbreite 128 ( ) dar. Die schwarz gepunktete Linie zeigt die
optimale Eﬃzienz.
3
Laufzeit in Sekunden
Laufzeit in Sekunden
3
2
1
0
2
1
0
1
20
40
60
80
100
120
140
1
20
40
Bandbreite
60
80
100
120
140
Bandbreite
Abbildung 3.8: Die Laufzeit der sequentiellen und parallelen SOFFT und iSOFFT für die Bandbegrenzung B = 128. Im linken Diagramm ist die Laufzeit der SOFFT zu sehen. Das rechte Diagramm zeigt
die Laufzeit der iSOFFT. Die rote Kurve stellt die Laufzeit der jeweiligen sequentiellen Version dar. Die
blaue Kurve entspricht der Laufzeit der jeweiligen parallelisierten Version. Die Anzahl der für die parallele SOFFT/iSOFFT verwendeten Einheiten beträgt 8.
eingetragen. Die schwarz gepunktete Linie stellt den optimalen Speedup eines parallelen Algorithmus dar. Die farbigen Kurven sind die Speedup-Kurven der SOFFT bzw. iSOFFT für
die verschiedenen Bandbreiten.
Abbildung 3.7 zeigt die Ergebnisse des Benchmarks zur Messung der Eﬃzienz aus Abschnitt
3.1. Die schwarz gepunktete Linie stellt die optimale Eﬃzienz eines parallelen Algorithmus
dar. Die farbigen Kurven entsprechen der Eﬃzienz der SOFFT bzw. iSOFFT für die verschiedenen Bandbreiten in Abhängigkeit von der Anzahl an Einheiten.
Ergebnisse der Benchmarks
B
16
32
64
128
256
512
1024
m = m0 = 0
m = B/2, m0 = 0
m = m0 = B/2
1.074883e − 14
2.201229e − 15
5.445416e − 15
3.304344e − 15
9.580270e − 15
2.375766e − 15
4.127737e − 14
7.338601e − 15
9.812105e − 15
8.944644e − 15
2.254168e − 15
2.981478e − 15
1.420379e − 13
1.677987e − 14
1.769254e − 14
2.176345e − 14
3.636223e − 15
3.849183e − 15
1.507377e − 13
4.969292e − 14
1.214394e − 13
1.632391e − 14
7.603809e − 15
1.861255e − 14
1.414312e − 12
1.367662e − 13
5.971587e − 13
1.081417e − 13
1.480205e − 14
6.473145e − 14
2.567475e − 12
3.908177e − 13
1.184365e − 12
1.389016e − 13
2.991182e − 14
9.075483e − 14
1.714953e − 11
1.100623e − 12
1.260088e − 08
6.564408e − 13
5.960425e − 14
6.827068e − 10
39
Tabelle 3.9: Rechengenauigkeit der DWT und iDWT, gemessen durch den in Abschnitt 3.1 beschriebenen Genauigkeits-Benchmark. Für jede Bandbreite gibt es zwei Zeilen, von denen die obere den
relativen Fehler und die untere den absoluten Fehler enthält.
In Abbildung 3.8 sind die Laufzeiten der sequentiellen SOFFT und iSOFFT im Vergleich
zur jeweiligen parallelisierten Version mit acht Einheiten zu sehen. Die roten Kurven stellen
dabei die Laufzeit der sequentiellen SOFFT/iSOFFT dar und die blauen Kurven die Laufzeit
der parallelisierten SOFFT/iSOFFT mit acht Einheiten.
Tabelle 3.9 zeigt abschließend die durch den Genauigkeits-Benchmark erzeugten Fehlerwerte für die DWT und iDWT. In der Tabelle sind die relativen und absoluten Fehler aufgeführt,
welche für die jeweiligen Zweierpotenzen der Bandbreite mittels Hin- und Rücktransformation durch DWT bzw. iDWT erzeugt wurden.
41
4
Diskussion und Schlussfolgerung
In diesem Kapitel wird auf die Messwerte aus Kapitel 3, Abschnitt 3.2 eingegangen, welche
durch die entsprechenden Benchmarks aus Abschnitt 3.1 erzeugt wurden. Insbesondere werden die Messwerte in Hinblick auf die Anwendbarkeit der parallelisierten SOFFT/iSOFFT in
der Praxis interpretiert. Im Anschluss daran wird ein Ausblick auf mögliche Weiterentwicklungen der vorgestellten Algorithmen gegeben.
Die farbigen Kurven in Abbildung 3.6 stellen die Speedup-Werte der parallelisierten SOFFT/
iSOFFT für die Bandbreiten 16, 32, 64 und 128 dar. Zusätzlich dazu ist der optimale Speedup schwarz gepunktet eingezeichnet. Ist ein Algorithmus perfekt parallelisiert, so beträgt
die Laufzeit dieses Algorithmus für p Einheiten Tp (n) = T ∗ (n)/p, wodurch der Speedup
dann
Sp (n) =
T ∗ (n)
T ∗ (n)
p
=
T ∗ (n) · p
=p
T ∗ (n)
beträgt. Angenommen, ein sequentieller Algorithmus benötigt zum Lösen eines Problems
mit einer Einheit eine Rechenzeit von T ∈ N. Dann erwartet man dementsprechend, dass
bei einer Nutzung von T Einheiten die Rechenzeit von T auf 1 sinkt, denn T Einheiten sind
in der Theorie T mal schneller als eine Einheit. Dies ergäbe also einen Speedup von T . In der
Praxis ist ein solcher Speedup jedoch nicht erreichbar. Es gibt zahlreiche Einflussfaktoren,
welche auf realen Systemen die perfekten Vorraussetzungen der Theorie stören. Bei Mehrprozessorsystemen müssen Einheiten z.B. miteinander kommunizieren. Der Master-Thread
wartet, wie in Kapitel 2, Abschnitt 2.2 erklärt, solange mit weiteren Berechnungen, bis alle
anderen Threads aus seinem Team wieder mit ihm vereinigt sind. Um zu wissen, ob noch andere Threads seines Teams arbeiten, muss er folglich mit diesen kommunizieren. Ein anderes
Problem besteht darin, dass eine Einheit nicht ausschließlich die Daten aus dem Arbeitsspeicher lädt, welche für die aktuellen Operationen benötigt werden. Um Rechenzeit zu sparen,
werden stattdessen ganze Speicherbereiche in den internen Speicher der Einheit, den sogenannten Cache geladen. Die Einheit ist in der Lage, auf diesen Speicher sehr viel schneller
zuzugreifen, als auf den Arbeitsspeicher. Arbeiten allerdings mehrere Einheiten mit denselben Daten und laden Speicherbereiche mit einer großen Überschneidung zu den Speicherbereichen der anderen Einheiten in den Cache, müssen während der Berechnung die Caches
synchronisiert werden. Der Grund hierfür liegt darin, dass Einheiten bei parallelen Berechnungen in gleichen Speicherbereichen eventuell mit nicht aktuellen Daten arbeiten. Diese
und andere Faktoren reduzieren den in der Praxis maximal erreichbaren Speedup.
Der gemessene Speedup in Abbildung 3.6 weist für eine geringe Anzahl an Einheiten ein
stärkeres Wachstum auf, als bei einer größeren Anzahl an verwendeten Einheiten. Es ist also
anzunehmen, dass der Speedup bei wachsender Zahl an Einheiten ab einem gewissen Punkt
42
Diskussion und Schlussfolgerung
nicht mehr steigt. Dies ist ein natürliches Phänomen, da der Overhead mit zunehmender
Anzahl an Einheiten zwangsläufig größer werden muss. Ein weiterer Grund ist, dass nicht
der ganze Algorithmus parallelisiert wurde. Aufgrund der verhältnismäßig kleinen Bandbreiten benötigt beispielsweise die parallele Ausführung der DWT eine längere Laufzeit als
die sequentielle Ausführung. Bei steigender Anzahl an Einheiten spielen ab einem gewissen Zeitpunkt ausschließlich die Laufzeit des sequentiellen Anteils der SOFFT/iSOFFT und
der Overhead eine Rolle. Der Overhead kann sich allerdings bereits früher bemerkbar machen. Für kleine Bandbreiten z.B. hat die sequentielle SOFFT/iSOFFT eine geringere Laufzeit als die parallele. Daher wurde in dieser Arbeit die SOFFT/iSOFFT erst ab der Bandbreite
B = 20 parallel ausgeführt.
Die Ergebnisse des Eﬃzienz-Benchmarks in Abbildung 3.7 bestätigen die Ergebnisse des
Speedup-Benchmarks. Mit Hilfe der Eﬃzienz kann eine Aussage darüber getroﬀen werden,
um wie viel schneller die Durchführung eines Algorithmus ist, wenn mehr als eine Einheit
verwendet wird. Eine Eﬃzienz von 1 für p Einheiten sagt aus, dass die Berechnungen mit p
Einheiten p mal schneller sind. Die optimale Eﬃzienz beträgt dementsprechend 1. In Abbildung 3.7 ist zu sehen, dass in dieser Arbeit beispielsweise für die Bandbreite B = 128 und
die Anzahl von p = 8 Einheiten eine Eﬃzienz von ca. 0.4 erreicht wurde. Dieser Wert impliziert, dass die parallele Ausführung der SOFFT/iSOFFT hier um etwa den Faktor 3 schneller
ist als die sequentielle. Dies lässt sich anhand von Abbildung 3.8 gut nachvollziehen. An dieser Stelle ist noch zu bemerken, dass in dieser Arbeit die DWT stets sequentiell durchgeführt
wurde, da sich dies trotz schlechterer Zeitkomplexität für alle verwendeten Bandbreiten als
vorteilhaft gegenüber der parallelen Ausführung erwiesen hat.
Zuletzt werden die Ergebnisse des Benchmarks für die Genauigkeit der DWT und iDWT
betrachtet. Die DWT bildet zusammen mit der FFT die Grundlage der SOFFT. Da die Rechengenauigkeit der FFT gut untersucht ist [Frigo und Johnson, 2012], ist der Benchmark
zur Messung des Transformations-Fehlers der DWT sehr wichtig. Vergleicht man die Werte in Tabelle 3.5 mit den entsprechenden Werten in [Kostelec und Rockmore, 2008] (dort
wurde ein sehr ähnlicher Test durchgeführt), so fällt auf, dass die DWT und iDWT dieser
Arbeit präziser arbeiten und dementsprechend geringere Fehler erzeugen. Selbst bei einer
Bandbreite von 256 ist der absolute Fehler nicht eindeutig von der Maschinenungenauigkeit
zu unterscheiden. Erst ab einer Bandbreite von 1024 entstehen größere Fehler. Diese Bandbreiten und größere Bandbreiten sind zum jetzigen Zeitpunkt allerdings nicht praxisrelevant,
weshalb dieser Verlust an Präzision von keiner großen Bedeutung ist.
Insgesamt zeigen die Ergebnisse der Benchmarks, dass die Parallelisierung die sequentielle
SOFFT und iSOFFT signifikant beschleunigt. Die in den vorherigen Kapiteln beschriebenen
Überlegungen und die daraus resultierende Parallelisierung der SOFFT und iSOFFT haben
somit einen praktischen Nutzen. Um weitere Verbesserungen der Laufzeit und somit des
43
Speedups und der Eﬃzienz zu erzielen, besteht wie in Kapitel 2 bereits erwähnt die Möglichkeit, die DWT und iDWT durch eine schnelle diskrete Polynom-Transformation zu ersetzen.
Durch diese Modifikation ergibt sich in der Theorie eine bessere Zeitkomplexität. Allerdings
steht hier die Untersuchung der Stabilität dieser modifizierten SOFFT und iSOFFT noch
aus. Eine weitere interessante Frage ist, inwieweit sich die in dieser Arbeit vorgestellte Parallelisierungsstrategie auch auf die SOFFT und iSOFFT mit nicht-äquidistanten Stützstellen
aus [Potts et al., 2009] übertragen lässt.
45
Abkürzungsverzeichnis
DSOFT diskrete Fourier-Transformation auf SO(3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
DWT diskrete Wigner-d-Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
EREW exclusive read exclusive write . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
FFT schnelle Fourier-Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
iDSOFT inverse diskrete Fourier-Transformation auf SO(3) . . . . . . . . . . . . . . . . . . . . . . . . . . 10
iDWT inverse diskrete Wigner-d-Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
iFFT inverse schnelle Fourier-Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
iSOFFT inverse schnelle Fourier-Transformation auf SO(3) . . . . . . . . . . . . . . . . . . . . . . . . . . 12
iSOFT inverse Fourier-Transformation auf SO(3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
MPI Message Passing Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
OpenMP Open Multi-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
ONB Orthonormalbasis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
O(3) orthogonale Gruppe im dreidimensionalen reellen Raum . . . . . . . . . . . . . . . . . . . . . . . . . . 4
SOFFT schnelle Fourier-Transformation auf SO(3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
SOFT Fourier-Transformation auf SO(3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
SO(3) spezielle orthogonale Gruppe im dreidimensionalen reellen Raum / Rotationsgruppe 4
47
Algorithmenverzeichnis
2.1
Sequentielle SOFFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
2.2
Sequentielle iSOFFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
2.7
Erster Umformungsschritt der SOFFT für die Parallelisierung . . . . . . . .
22
2.9
Zweiter Umformungsschritt der SOFFT für die Parallelisierung . . . . . . .
24
2.10 Parallelisierte SOFFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
2.11 Parallelisierte iSOFFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
49
Abbildungsverzeichnis
2.3 Beispiel einer Parallelisierung dargestellt als Fork-Join-Modell . . . . . . . .
17
2.4 Graphische Darstellung der Forster-Methodik . . . . . . . . . . . . . . . . .
19
2.5 Parallelisierung der äußeren For-Schleife . . . . . . . . . . . . . . . . . . .
20
2.6 Parallelisierung der inneren For-Schleife . . . . . . . . . . . . . . . . . . . .
20
2.8 Geometrischer Lösungsansatz der For-Schleifen-Optimierung . . . . . . . .
23
3.1 Benchmark zu Messung des Speedups (SOFFT) . . . . . . . . . . . . . . . .
34
3.2 Benchmark zur Messung des Speedups (iSOFFT) . . . . . . . . . . . . . . .
34
3.3 Benchmark zur Messung der Eﬃzienz (SOFFT) . . . . . . . . . . . . . . . .
35
3.4 Benchmark zur Messung der Eﬃzienz (iSOFFT) . . . . . . . . . . . . . . .
35
3.5 Benchmark zur Messung der Rechengenauigkeit der DWT aus Satz 1.17. . .
36
3.6 Ergebnisse des Speedup Benchmarks . . . . . . . . . . . . . . . . . . . . . .
37
3.7 Ergebnisse des Eﬃzienz-Benchmarks . . . . . . . . . . . . . . . . . . . . .
38
3.8 Laufzeit der sequentiellen/parallelisierten SOFFT/iSOFFT . . . . . . . . . .
38
3.9 Rechengenauigkeit der DWT . . . . . . . . . . . . . . . . . . . . . . . . . .
39
51
Literaturverzeichnis
Bülow, H. und Birk, A. (2011a). Spectral registration of noisy sonar data for underwater 3D
mapping. Auton. Robots, 30(3):307–331.
Bülow, H. und Birk, A. (2011b). Spectral registration of volume data for 6-DOF spatial transformations plus scale. In Proc. IEEE ICRA 2011, pages 3078–3083, Shanghai, China.
Cormen, T. H., Leiserson, C. E., Rivest, R. L., und Stein, C. (2009). Introduction to Algorithms,
Third Edition. The MIT Press, Cambridge, MA, USA.
Driscoll, J. R., Healy, D. M., und Rockmore, D. N. (1997). Fast Discrete Polynomial Transforms with Applications to Data Analysis for Distance Transitive Graphs. SIAM J. Comput.,
26(4):1066–1099.
Edmonds, A. R. (1996). Angular Momentum in Quantum Mechanics. Princeton University
Press, Princeton, NJ, USA.
Fischer, G. (2008). Lehrbuch der Algebra: mit lebendigen Beispielen, ausführlichen Erläuterungen und zahlreichen Bildern. Vieweg Verlag, Wiesbaden.
Frigo, M. und Johnson, S. G. (2012).
http://www.fftw.org/fftw3.pdf.
[Online; Zugriﬀ am
12. September 2015].
JáJá, J. (1992). An Introduction to Parallel Algorithms. Addison Wesley Longman Publishing
Co., Inc., Redwood City, CA, USA.
Koecher, M. (1997). Lineare Algebra und analytische Geometrie. Springer-Verlag Berlin Heidelberg, Berlin, Heidelberg.
Kostelec, P. J. und Rockmore, D. N. (2008). FFTs on the Rotation Group. J. Fourier Anal.
Appl., 14(2):145–179.
Kovacs, J. A. und Wriggers, W. (2002). Fast rotational matching. Acta Cryst., D58:1282–1286.
Matsumoto, M. und Nishimura, T. (1998). Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator. ACM Trans. Model. Comput.
Simul., 8(1):3–30.
Modersitzki, J. (2009). FAIR: Flexible Algorithms for Image Registration. Society for Industrial
and Applied Mathematics (SIAM), Philadelphia, PA, USA.
Potts, D., Prestin, J., und Vollrath, A. (2009). A fast algorithm for nonequispaced Fourier
transforms on the rotation group. Numer. Algorithms, 52(3):355–384.
Quinn, M. J. (2003). Parallel Programming in C with MPI and OpenMP. The McGraw-Hill
Companies, Inc., New York, NY, USA.
52
Literaturverzeichnis
Ritchie, D. W. (2011). High Performance Algorithms for Molecular Shape Recognition. Scientific Memoir, Université Henri Poincaré.
Roosta, S. H. (2000). Parallel Processing and Parallel Algorithms. Springer-Verlag New York,
Inc., New York, NY, USA.
Slabaugh, G., Fang, T., McBagonluri, F., Zouhar, A., Melkisetoglu, R., Xie, H., und Ünal,
G. (2008).
3-D shape modeling for hearing aid design.
IEEE Signal Process. Mag.,
25(5):98–102.
The GCC-Team (2015).
Optimize-Options.
https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#
[Online, Zugriﬀ am 01. Oktober 2015].
The Intel Corporation (2015).
http://ark.intel.com/de/products/75122/Intel-Core-i7-
4770-Processor-8M-Cache-up-to-3_90-GHz.
[Online; Zugriﬀ am 12. September 2015].
The OpenMP Architecture Review Board (2013).
OpenMP4.0.0.pdf.
http://www.openmp.org/mp-documents/
[Online; Zugriﬀ am 12. September 2015].

Download Report