Das Programm NKOERPER löst die Bewegungsgleichungen von Teilchen, die sich
unter dem Einfluß der wechselseitigen Gravitationskräfte bewegen. Es
benutzt dazu ein Runge-Kutta-Verfahren mit Adaption der Schrittweite. Die
entsprechenden Routinen runge_kutta und integriere sind Adaptionen von
Programmen aus dem Buch ''Theoretische Physik mit dem Personal Computer''
von Schmid, Spitz und Lösch.
Das Hauptprogramm main ist folgendermaßen aufgebaut:
Mit io_init werden alle Eingabe-Parameter wie Zahl der Teilchen,
ihre
Massen, Anfangsorte und -geschwindigkeiten eingelesen und über den
globalen IO-Block zurückgegeben.19 Die Routine dvector besorgt den Speicherplatz für Orts- und
Geschwindigkeits-Vektoren, die dann auf die Anfangsbedingungen gesetzt
werden. Nach diesen Vorbereitungen startet die Rechnung: Mit integriere
werden Ort und Geschwindigkeit der Teilchen für den gegenwärtigen
Zeitpunkt bestimmt und mit write_vector ausgegeben.
set_vector setzt die
Vektoren für die nächste Iteration und die Zeit wird um einen Schritt
erhöht. Nach der Rechnung werden noch mit io_close
abschließende Arbeiten
erledigt (Files geschloßen).
''integriere'' implementiert die Schrittweiten-Anpassung und ruft für die
eigentliche Integration die Routine runge_kutta auf, die
allgemeine Systeme
von Differentialgleichungen erster Ordnung lösen kann. Das hier benötigte
spezielle System
[xdot]a = [v]a (a: Teilchenindex)
[vdot]a = sum (b!=a) [m]b /([r]ab ** 3) * ([x]b - [x]a)
wird durch ''kepler_dgl'' definiert.
- Version 2.1
-
- Änderungen:
-
- Meßergebnisse:
-
time(-no): |
438 |
. | 2s |
time(-O0): |
364 |
. | 9s |
time(-O1): |
314 |
. | 7s |
time(-O2): |
331 |
. | 7s |
- CXpa(-O1):
-
|
CPU Time |
. | CPU Time |
Times |
. | Routine |
|
|
|
(less children) |
. | (plus children) |
Called |
. | Name |
|
|
|
306 |
. | 987 |
56 |
. | 0% |
546 |
. | 519 |
99 |
. | 7% |
3924 |
kepler_dgl |
|
239 |
. | 526 |
43 |
. | 7% |
539 |
. | 526 |
43 |
. | 7% |
39632727 |
mthd_pown |
|
1 |
. | 307 |
0 |
. | 2% |
547 |
. | 845 |
99 |
. | 9% |
981 |
runge_kutta |
- Analyse:
-
- NKOERPER vektorisiert nicht
- Zeit geht in kepler_dgl und sein ''Kind'' mthd_pown
- Bemerkungen:
-
- mthd_pown ist die Potenz-Funktion aus der Standard-Bibliothek. Da
nur die Potenz 3/2 auftritt, kann sie durch die Wurzel-Funktion ersetzt
werden, die auf der C3 in Hardware realisiert ist.
- Version 2.2
-
- Änderungen:
-
- pow-Funktion in kepler_dgl durch sqrt ersetzt
- Meßergebnisse:
-
time(-O1): |
168 |
. | 2s |
time(-O2): |
179 |
. | 3s |
- CXpa(-O1):
-
|
CPU Time |
. | CPU Time |
Times |
. | Routine |
|
|
|
(less children) |
. | (plus children) |
Called |
. | Name |
|
|
|
166 |
. | 538 |
99 |
. | 0% |
166 |
. | 546 |
99 |
. | 0% |
3924 |
kepler_dgl |
- Analyse:
-
- großer Zeitanteil der Potenz-Routine völlig weggefallen
- kepler_dgl vektorisiert nicht
- Bemerkungen:
-
- Die Schleifen in kepler_dgl sind zu groß und zu kompliziert, sie
werden in Teile aufgespalten.
- Man kann die Symmetrie der Kräfte (F_ab =
F_ba, actio = reactio) mit
folgendem Trick ausnutzen: Statt des Vektors rinv3[b] in Z.92, der nur
für einen Wert von a (ein Teilchen) gilt, wird ein Array rinv3[a][b]
berechnet, wobei dann die Symmetrie eingehen kann.
- Zur einfachen dynamischen Erzeugung eines Arrays wird eine Funktion
darray analog zu dvector benötigt. Damit auf das Array mit der
üblichen Index-Notation a[i][j] zugegriffen werden kann, muß man einen
Vektor a von Zeigern auf die Zeilen des Arrays erzeugen und entsprechend
initialisieren.
- Version 2.3
-
- Änderungen:
-
- Hilfsarray rinv3[a][b] eingeführt
- Symmetrie ausgenutzt durch Umstrukturierung der Schleifen
- neue Routine darray in vektor_utl.c
- Meßergebnisse:
-
time(-O1): |
128 |
. | 3s |
time(-O2): |
142 |
. | 2s |
- CXpa:
-
|
CPU Time |
. | CPU Time |
Times |
. | Routine |
|
|
|
(less children) |
. | (plus children) |
Called |
. | Name |
|
|
|
126 |
. | 037 |
98 |
. | 7% |
126 |
. | 082 |
98 |
. | 7% |
3924 |
kepler_dgl (-O1) |
|
140 |
. | 455 |
98 |
. | 9% |
140 |
. | 463 |
98 |
. | 9% |
3924 |
kepler_dgl (-O2) |
- Optimierungsbericht für kepler_dgl:
-
Line |
Iter. |
Reordering |
Optimizing / Special |
Exec. |
Num. |
Var. |
Transformation |
Transformation |
Mode |
7 |
a |
Scalar |
|
|
72 |
b |
Scalar |
|
|
77 |
i |
FULL VECTOR |
Reduction |
|
|
|
|
|
|
87 |
a |
Scalar |
|
|
89 |
b |
Scalar |
|
|
|
|
|
|
|
96 |
a |
Scalar |
|
|
98 |
i |
Scalar |
|
|
103 |
b |
Scalar |
|
|
|
|
|
|
|
Line |
Iter. |
Analysis |
|
|
Num. |
Var. |
|
|
|
70 |
a |
Inner loop has varying trip count |
|
|
72 |
b |
Unable to distribute |
|
|
87 |
a |
Inner loop has varying trip count |
|
|
89 |
b |
Insufficient vector code |
|
|
96 |
a |
Unable to distribute |
|
|
98 |
i |
Unable to distribute |
|
|
103 |
b |
Insufficient vector code |
|
|
- Analyse:
-
- Skalarzeit wesentlich verbessert
- einzige vektorisierte Schleife: Z.77
- Bemerkungen:
-
- Schleife in Z.77 hat nur Vektorlänge d=2, sollte also besser nicht
vektorisieren. Da d aber erst zur Laufzeit bestimmt wird, kann dies kein
Compiler erkennen.
- Application-Compiler bringt gar nichts.
- kepler_dgl enthält zwei Teile: ``x' = v``
und ``v' = F``. Der erste Teil
kann in einer eigenen Schleife über alle Teilchen bestimmt werden.
- Um die Extrabehandlung der Diagonalen (a==b) in Z.105 - Z.108 mit dem
continue loszuwerden, setzen wir bei der Berechnung von rinv3
rinv3[a][a] = 0.
- Version 2.4
-
- Änderungen:
-
- Schleifenstruktur in kepler_dgl vereinfacht
- Meßergebnisse:
-
time(-O1): |
122 |
. | 2s |
time(-O2): |
120 |
. | 9s |
- CXpa:
-
|
CPU Time |
. | CPU Time |
Times |
. | Routine |
|
|
|
(less children) |
. | (plus children) |
Called |
. | Name |
|
|
|
120 |
. | 066 |
98 |
. | 7% |
120 |
. | 074 |
98 |
. | 7% |
3924 |
kepler_dgl (-O1) |
|
119 |
. | 051 |
98 |
. | 7% |
119 |
. | 059 |
98 |
. | 7% |
3924 |
kepler_dgl (-O2) |
- Analyse:
-
- keine Vektorisierung auch einfacher Schleifen in
kepler_dgl, da das
dynamische Array rinv3 mit Zeilenpointern aufgebaut ist.
- Bemerkungen:
-
- Der Zugriff auf Arrays über Zeilenzeiger verhindert, daß Spalten
als Vektoren erkannt werden können. Um Vektorisierung überhaupt
möglich zu machen, müssen die dynamischen Arrays also anders erzeugt
werden. Dazu wird ein Array als eine Datenstruktur mit einem langen
Vektor und Dimensionen definiert und über ein spezielles INDEX-Makro
darauf zugegriffen. Die Elemente werden dann statt mit der Notation
a[i][j] durch INDEX(a,i,j) angesprochen.
- Version 2.5
-
- Änderungen:
-
- Typ dmatrix und Makro INDEX in nkoerper.h definiert
- darray in vektor_utl.c angepaßt
- kepler_dgl angepaßt
- Meßergebnisse:
-
time(-O1): |
101 |
. | 6s |
time(-O2): |
107 |
. | 7s |
- Optimierungsbericht für kepler_dgl:
-
Line |
Iter. |
Reordering |
Optimizing / Special |
Exec. |
Num. |
Var. |
Transformation |
Transformation |
Mode |
71 |
a |
Scalar |
|
|
73 |
b |
Scalar |
|
|
77 |
i |
FULL VECTOR |
Reduction |
|
|
|
|
|
|
87 |
a |
Scalar |
|
|
89 |
b |
Scalar |
|
|
|
|
|
|
|
96 |
a |
FULL VECTOR |
|
|
|
|
|
|
|
103 |
i |
Scalar |
|
|
|
|
|
|
|
109 |
a |
Scalar |
|
|
111 |
i |
Scalar |
|
|
114 |
b |
Scalar |
|
|
|
|
|
|
|
Line |
Iter. |
Analysis |
|
|
Num. |
Var. |
|
|
|
71 |
a |
Inner loop has varying trip count |
|
|
73 |
b |
Unable to distribute |
|
|
87 |
a |
Inner loop has varying trip count |
|
|
89 |
b |
Insufficient vector code |
|
|
103 |
i |
Insufficient vector code |
|
|
109 |
a |
Unable to distribute |
|
|
111 |
i |
Unable to distribute |
|
|
114 |
b |
Insufficient vector code |
|
|
Loop Performance Analysis |
|
For kepler_dgl.c:kepler_dgl |
|
|
|
Line |
Times |
Iteration |
Count |
|
Number |
Exec. |
Min |
Max |
Avg |
Total CPU Time |
|
71 |
3924 |
101 |
101 |
101 |
109 |
. | 110868 |
73 |
392400 |
1 |
100 |
50 |
108 |
. | 079641 |
77 |
19816200 |
2 |
2 |
2 |
44 |
. | 060601 |
87 |
3924 |
101 |
101 |
101 |
10 |
. | 615843 |
89 |
392400 |
1 |
100 |
50 |
9 |
. | 681291 |
96 |
3924 |
101 |
101 |
101 |
0 |
. | 019290 |
103 |
3924 |
202 |
202 |
202 |
0 |
. | 223310 |
109 |
3924 |
101 |
101 |
101 |
61 |
. | 225854 |
111 |
396324 |
2 |
2 |
2 |
60 |
. | 099856 |
114 |
792648 |
101 |
101 |
101 |
57 |
. | 622604 |
- Analyse:
-
- Schleife Z.96 (rinv3[a][a] = 0) vektorisiert, d.h. grundsätzlich ist
Vektorisierung mit der neuen Datenstruktur möglich
- übrige Schleifen enthalten alias-Probleme
- Hauptschleifen: Z.71-84 und Z.114
- Bemerkungen:
-
- Die Schleifen Z.71-84 können durch Einführung eines Hilfsarrays
r[a][b], das die Abstände der Teilchen enthält, vereinfacht werden.
- Die Schleifen über Teilchenindizes müssen explizit nach innen
gebracht werden, die über Raumindizes nach außen, weil der Compiler
nicht wissen kann, daß die Teilchenzahl wesentlich größer ist als die
Dimension des Raums.
- Version 2.6
-
- Änderungen:
-
- Vereinfachung und Aufteilung der Schleifenstruktur in Z.71-84
- Meßergebnisse:
-
time(-O1): |
100 |
. | 5s |
time(-O2): |
88 |
. | 4s |
- Optimierungsbericht für kepler_dgl:
-
Line |
Iter. |
Reordering |
Optimizing / Special |
Exec. |
Num. |
Var. |
Transformation |
Transformation |
Mode |
74 |
a |
Scalar |
|
|
76 |
b |
FULL VECTOR |
|
|
|
|
|
|
|
82 |
a |
Scalar |
|
|
84 |
i |
Scalar |
|
|
86 |
b |
Scalar |
|
|
|
|
|
|
|
93 |
a |
Scalar |
|
|
95 |
b |
Scalar |
|
|
|
|
|
|
|
103 |
a |
Scalar |
|
|
105 |
b |
Scalar |
|
|
|
|
|
|
|
112 |
a |
FULL VECTOR |
|
|
|
|
|
|
|
120 |
i |
Scalar |
|
|
|
|
|
|
|
126 |
a |
Scalar |
|
|
128 |
i |
Scalar |
|
|
131 |
b |
Scalar |
|
|
- Analyse:
-
- kaum Vektorisierung wegen Alias-Problem (z.B. zwischen r und x)
- Bemerkungen:
-
- Auch der Application-Compiler schafft es nicht, die Aliasse
aufzulösen. Anscheinend versagt die Pointer-Analyse bei dynamisch
erzeugten Arrays.
- Version 2.7
-
- Änderungen:
-
- no_recurrence-Direktive eingefügt in Z.86, 96, 107, 123
- Meßergebnisse:
-
- Optimierungsbericht für kepler_dgl:
-
Line |
Iter. |
Reordering |
Optimizing / Special |
Exec. |
Num. |
Var. |
Transformation |
Transformation |
Mode |
74 |
a |
Scalar |
|
|
76 |
b |
FULL VECTOR |
|
|
|
|
|
|
|
82 |
a |
Scalar |
|
|
84 |
i |
Scalar |
|
|
87 |
b |
FULL VECTOR |
|
|
|
|
|
|
|
94 |
a |
Scalar |
|
|
97 |
b |
FULL VECTOR |
|
|
|
|
|
|
|
105 |
a |
Scalar |
|
|
108 |
b |
FULL VECTOR |
|
|
|
|
|
|
|
115 |
a |
FULL VECTOR |
|
|
|
|
|
|
|
124 |
i |
FULL VECTOR |
|
|
|
|
|
|
|
130 |
a |
Scalar |
|
|
132 |
i |
Scalar |
|
|
135 |
b |
Scalar |
|
|
- Analyse:
-
- Schleifen in Z.74-124 voll vektorisiert
- Bemerkungen:
-
- Nun müssen noch die Schleifen Z.130-135 zerlegt werden.
- Version 2.8
-
- Änderungen:
-
- extra Schleife für fv=0 eingeführt
- Schleife über Teilchenindex a nach innen getauscht
- no_recurrence-Direktive vor Schleife über a
- Meßergebnisse:
-
time(-O1): |
93 |
. | 8s |
time(-O2): |
22 |
. | 3s |
- CXpa:
-
|
CPU Time |
. | CPU Time |
Times |
. | Routine |
|
|
|
(less children) |
. | (plus children) |
Called |
. | Name |
|
|
|
20 |
. | 708 |
93 |
. | 0% |
20 |
. | 723 |
93 |
. | 0% |
3924 |
kepler_dgl |
- Analyse:
-
- volle Vektorisierung aller möglichen Schleifen in
kepler_dgl
- Vektorisierungs-Beschleunigung 4.2
- Bemerkungen:
-
- Da kepler_dgl auch bei guter Vektorisierung immer noch den größten
Teil der Zeit verbraucht, bleibt zur weiteren Beschleunigung nur die
Möglichkeit, den Algorithmus zum Lösen der DGL mit variabler
Schrittweite so zu verbessern, daß bei gleicher Genauigkeit weniger oft
kepler_dgl aufgerufen wird. Dazu wollen wir uns die entsprechenden
NAG-Routinen vornehmen.
- Version 3.1
-
- Änderungen:
-
- integriere nur noch Interface zur NAG-Routine D02BAF
- in kepler_dgl die Variablen x und v zu einem Vektor (x,v) zusammengefaßt
- Meßergebnisse:
-
- CXpa:
-
|
CPU Time |
. | CPU Time |
Times |
. | Routine |
|
|
|
(less children) |
. | (plus children) |
Called |
. | Name |
|
|
|
2 |
. | 435 |
92 |
. | 7% |
2 |
. | 437 |
92 |
. | 8% |
495 |
kepler_dgl |
- Analyse:
-
- kepler_dgl wird statt 3924 mal nur 495 mal aufgerufen
- Bemerkungen:
-
- Um das Ergebnis beurteilen zu können, müssen auch die Genauigkeiten
beider Verfahren verglichen werden. Die Konstante eps hat laut NAG-Manual
eine ähnliche Bedeutung wie in Version 2.1 - 2.8, allerdings ist nicht
klar, was dies für die Genauigkeit des Ergebnisses bedeutet. Im Rahmen
der Ausgabe-Genauigkeit (4 Dezimalen) sind beide Verfahren identisch, so
daß sich die Frage stellt, welche Genauigkeit überhaupt nötig ist und
ob man nicht auch dadurch enorm CPU-Zeit sparen kann, daß man die
Genauigkeit heruntersetzt. Dies erfordert aber weitergehende
Untersuchungen, die hier nicht angestellt werden.
- Dem NAG-Manual kann man entnehmen, daß die Routine D02BAF im
Normalfall, insbesondere bei nicht zu zeitaufwendiger rechter Seite
anzuwenden ist. Ansonsten wird D02CAF (''Adams-Verfahren'') empfohlen. Als
Änderung ist nur ein Austausch des Routinennamnes und der Größe des
Workspace nötig.
- Version 3.2
-
- Änderungen:
-
- NAG-Routine D02BAF durch D02CAF ersetzt
- Meßergebnisse:
-
- CXpa:
-
|
CPU Time |
. | CPU Time |
Times |
. | Routine |
|
|
|
(less children) |
. | (plus children) |
Called |
. | Name |
|
|
|
0 |
. | 787 |
79 |
. | 8% |
0 |
. | 788 |
79 |
. | 9% |
142 |
kepler_dgl |
- Analyse:
-
- Zeitgewinn durch drastisch reduzierte Zahl der Aufrufe von
kepler_dgl
- Bemerkungen:
-
- Beschleunigung gegenüber Version 2.1: 350 !
- Zusammenfaßung der Optimierungen:
-
Version |
time(-O2)/s |
Änderungen |
|
2.1 |
331 |
. | 7 |
Ausgangsversion |
2.2 |
179 |
. | 3 |
pow-Funktion in kepler_dgl durch sqrt ersetzt |
2.3 |
142 |
. | 2 |
Symmetrie, Vereinfachung der Schleifen in kepler_dgl |
2.4 |
120 |
. | 9 |
Schleifenstruktur in kepler_dgl weiter vereinfacht |
2.5 |
107 |
. | 7 |
Typ dmatrix und Makro INDEX eingeführt |
2.6 |
88 |
. | 4 |
Vereinfachung der Schleifenstruktur in Z.71-84 |
2.7 |
65 |
. | 4 |
no_recurrence-Direktiven |
2.8 |
22 |
. | 3 |
Vereinfachung der Schleifen in Z.130-135, Direktive |
3.1 |
2 |
. | 7 |
integriere durch NAG-Routine D02BAF ersetzt |
3.2 |
0 |
. | 9 |
NAG-Routine D02CAF statt D02BAF |
Peter Junglas 18.10.1993