-
Die Erfindung betrifft ein Verfahren
zum Betrieb eines Vektorrechners, um eine Gleichung vom Typ: A q = r (1)zu lösen, wobei
A eine m × m
Matrix ist und q und r jenweils Vektoren der Dimension m sind.
-
Insbesondere betrifft die Erfindung
ein Verfahren zur Lösung
der Gleichungen (Δ – L) p =
r (2)und Δ–1 (Δ – U) q =
p (3)wobei
L eine untere m × m
Matrix, U eine obere m × m
Matrix und Δ eine
m × m
Diagonal-Matrix ist, und wobei p, q und r jeweils Vektoren der Dimension
m sind.
-
Gleichungen des obengenannten Typs
treten bei Differentialgleichungen bei der Lösung von Randwertproblemen,
insbesondere bei elliptischen partiellen Differentialgleichungen,
auf. Diese Randwertprobleme ergeben sich bei einer Vielzahl numerischer
Berechnungen technischer Größen, so
etwa bei der Berechnung der Temperaturverteilung in einem dreidimensionalen
Körper
oder bei der Berechnung einer Fluidströmung um eine gegebene geometrische
Form. Ein weiteres einfaches Beispiel ist die Lösung der Poisson-Gleichung
mit Dirichlet-Randbedingungen.
-
Üblicherweise
wird für
die numerische Lösung
das Problem diskretisiert, indem der Raum, in dem die Lösung zu
berechnen ist, mit einem Punktgitter überzogen wird, und wobei dann
numerisch eine Lösung
für jeden
Gitterpunkt berechnet wird. Hierzu kann die Lösungsfunktion um jeden Gitterpunkt
beispielsweise in eine Taylorreihe entwickelt werden, und dann können die
Differentialquotienten durch Differenzquotienten ersetzt werden.
Zur Berechnung einer Näherungslösung an
einem Gitterpunkt reicht dann die Kenntnis der Näherungslösung an den benachbarten Gitterpunkten
aus, wenn die Restglieder der Reihenentwicklung vernachlässigt werden
können.
Ein solches Verfahren ist beispielsweise in "Numerische Mathematik für Ingenieure
und Physiker" W.
Törnig,
Band 2, Springer-Verlag Berlin offenbart.
-
Beispielsweise kann es in einem zweidimensionalen
Raum mit einem quadratischen Gitter ausreichen, jeweils die Lösung an
den unmittelbar benachbarten vier Gitterpunkten zu kennen, um die
Näherungslösung an
dem von diesen vier Gitterpunkten umgebenen Gitterpunkt zu berechnen.
-
1 zeigt
diesen Fall, indem ein sog. 5-Punkt-Stencil verwendet wird, während 2 den entsprechenden Fall
für den
dreidimensionalen Raum zeigt, wobei sich dann ein 7-Punkt-Stencil ergibt.
-
Mit dem oben beschriebenen Ansatz
läßt sich
das Problem bei geeigneter Numerierung der Gitterpunkte für den gezeigten
7-Punkt-Stencil wie folgt darstellen:
aijk,1 qijk-1 +
aijk,2 qij-1k +
aijk,3 qi-1jk +
aijk,4 qijk + aijk,5 qi+1jk + aijk,6 qij+1k + aijk,7 qijk+1 = rijk
(4)
-
Für
jeden Gitterpunkt (ijk) aus der Menge G mit G = {(ijk) 1 ≤ i ≤ nx, 1 ≤ j ≤ ny, 1 ≤ k ≤ nz} (5)wobei
nx, ny und nz jeweils die Anzahl der Gitterpunkte in
x, y bzw. z-Richtung sind.
-
Die Matrix A enthält bei geeigneter Wahl der Numerierung
der Gitterpunkte für
gewöhnlich
nur wenige Außerdiagonalelemente,
die von Null verschieden sind.
-
Die obengezeigte Gleichung 1 wird
dann nicht direkt gelöst,
sondern zunächst
unter Ausnutzung der Beziehung A = D – L – U (6)und eines
Präkonditionierers M = P1 × P2
(7)umgeschrieben
in (P1
–1 A
P2
–1) (P2 q)
= p1
–1 r (8)
-
Hierbei ist D die Diagonalmatrix
aus den Elementen der Diagonalen der ursprünglichen Matrix A.
-
Ein oft verwendeter und erfolgreicher
Präkonditionierer
ist der sog. "Incomplete
lower upper preconditioner" (ILU-Preconditioner). M = (∆ – L') ∆–1 (∆ – U') (9)
-
Hierbei ist ∆ = D – Diagonalelemente (L Δ–1 U) – α Zeilensumme
(Außer-Diagonalelemente
(L Δ–1 U))
mit 0 ≤ α ≤ 1; Außer-Diagonalelemente
(F) = F – Diagonalemente
(F); und Zeilensumme bedeutet eine Diagonalmatrix, bei der jedes
Element der Diagonalen der Summe der Elemente derselben Zeile entspricht.
Dies ist die bekannte Gustafsson-Beschleunigung (I. Gustafsson "A Class of First
Order Factorization Methods",
in BIT, ISSN 0006-3835, Vol. 18, No. 2 978), Seiten 142-156). Die
Matrizen U' und L' sind in Gleichung
(9) mit einem Apostroph gekennzeichnet, um darauf hizuweisen, daß sie nicht
notwendigerweise identisch zu den Matrizen U und L aus Gleichung
(6) sind, obwohl dies beim beschriebenen ILU-Ansatz der Fall ist. Geeignete Änderungen
der Matrizen U' und
L'können zu
einer verbesserten Konvergenz führen.
-
Zur Bestimmung von Δ in Gleichung
(9) gibt es weitere Ansätze,
etwa den SSOR-Ansatz mit Δ = D.
-
Die Präkonditionierung entspricht
der Lösung
des linearen Gleichungssystems M q = r (10)wobei
r der bekannte Vektor ist. Eine schematische Ansicht dieser Gleichung
ist in 3A gezeigt.
-
Die Lösung kann durch Lösen der
Gleichungen (Δ – L) p =
r (11)und Δ–1 (Δ – U) q =
p (12)erhalten
werden.
-
Wenn das anfängliche Problem auf ein strukturiertes
Gitter gerichtet war, können
die von Null verschiedenen Elemente der Matrizen U und L regelmäßig angeordnet
werden. Zum Beispiel sind im Fall des 7-Stencils im regelmäßigen dreidimensionalen Gitter
mit jeweils nx, ny und
nz Gitterpunkten und lexikographischer Numerierung
der Gitterpunkte in der Matrix L nur in der ersten Nebendiagonalen,
in der nx-ten Nebendiagonalen und in der
(nx ny)-ten Nebendiagonalen
von Null verschiedene Elemente vorhanden. In anderen Worten sind
nur die Elemente an Li,i-1,
Li,i-nx und Li,i-nx ny von
Null verschieden.
-
Die Gleichungen (11) und (12) können durch entsprechende
Vorwärts-
bzw. Rückwärtssubstitutionen
gelöst
werden. p = ∆–1 (r
+ L p) (13)
q = P + ∆–1 U
q (14)
-
Ein Algorithmus (1) für diese
Substitutionen ist: für
i = 1, ..., m
pi = (1/Δi,i)
(ri + Σj<iLi,j pj)
Ende
der Schleife für
i
für
i = m, ..., 1
qi = pi +
(1/Δi,i) Σj>iUi,j qj
Ende
der Schleife für
i.
-
Dieses Verfahren ist als Hyper-Plane-Verfahren
bekannt. Man kann zeigen, daß Hyperebenen (Hyper-Planes)
aus Gitterpunkten konstruiert werden können, die unabhängig voneinander
berechnet werden können,
da die Lösung
an einem beliebigen Gitterpunkt der Hyperebene nicht von der Lösung irgendeines
anderen Punktes derselben Ebene abhängt. Entsprechende Hyperebenen
für den
dreidimensionalen Raum zeigt 3B.
Für den 7-Punkt-Stencil
gilt:
H(l) = {i =
ix + (iy – 1)nx + (iz – 1) nx ny | ix +
iy + iz = l + 2} (15)
-
Der Nachteil dieses Verfahrens besteht hauptsächlich darin,
daß die
Effizienz einer Warteschlangen-Verarbeitung bzw. einer vektoriellen
Verarbeitung in einem Vektorprozessor durch die indirekte Adressierung
in dem oben beschriebenen Algorithmus (1) stark vermindert wird.
Vorteilhafter wäre es,
immer mit einem konstanten Speicherabstand (stride) auf den Speicher
zugreifen zu können
(siehe auch "Computer
Architecture"; Computational
Science Education Project, Abschnitt 3.3.1 "Interleaved Memory"). Eine Methode zur Erhöhung der
Speichergeschwindigkeit bei Vektorzugriffen durch lineare Adresstransformationen
wird von D. T. Harper III in „Increased
Memory Performance During Vector Accesses Through the Use of Linear
Address Transformations" in
IEEE Transactions on Computers, ISSN 0018-9340, Vol. 41, No. 2,
S. 227-230 (1992) beschrieben.
EP 0 479 235 A2 offenbart einen Vektorrechner,
der in einer indirekten oder einer Stride-Vektorrechenbetriebsart
arbeiten kann.
-
Außerdem ist das Verfahren nicht
für Rechner
mit geteilten parallelen Speicherblöcken und mehreren Vektorprozessoren
geeignet, da keine Parallelität
bezüglich
der Größe 1 der
Hyperebene auftritt, und da viele Hyperebenen H(l) alleine zu klein sind,
um sinnvoll auf mehrere Prozessoren verteilt zu werden. Des weiteren
ist bei den genannten Hyper-Plane-Verfahren keine Verallgemeinerung
auf höherdimensionale
Räume oder
andere als die bekannten 5- bzw. 7-Punkt-Stencils bekannt. Ein Ansatz,
diesen Problemen zu begegnen, ist es, die Gleichungen 11 und 12
nicht exakt zu lösen,
sondern durch Iteration eine Näherungslösung zu
berechnen. Üblicherweise
wird hierfür
eine von-Neumann-Reihenentwicklung verwendet. p = (Δ – L)–1 r
= Δ–1 (I – L Δ–1)–1 r
= Δ–1 r
+ Δ–1 L Δ–1 r
+ (Δ–1 L)2 Δ–1 r
+ (Δ–1 L)3 Δ–1 r
+ ... + (Δ–1 L)ntr Δ–1 r (16).
-
Dieses Verfahren ist dann zwar voll
vektorisierbar und kann parellel bearbeitet werden, da immer nur
Matrix-Vektor-Produkte berechnet werden, um jedoch eine Näherungslösung zu
erhalten, die in etwa so gut wie jene des Hyper-Plane-Verfahrens ist, muß meistens
die von-Neumann-Reihe bis zu den Gliedern dritter oder vierter Ordnung
berechnet werden. Der Rechenaufwand ist dann 3 bis 4 mal so groß, was dazu führen kann,
daß der
Vorteil aus der Vektorisierung und der Parallelisierung aufgehoben wird.
-
Ein weiterer Ansatz, um das Hyper-Plane-Verfahren
in geeigneter Weise für
Vektorrechner abzuwandeln, ist von Takumi Washio und Ken Hayami
in "Overlapped Multicolor
MILU Preconditioning", in
SIAM Journal on Scientific Computing, ISSN 1064-8275, Band 16, Nr.
(1995), Seiten 636-651 offenbart worden.
-
Entsprechend diesem Ansatz werden
für den Fall
des dreidimensionalen Raumes und des 7-Punkt-Stencils zunächst die
Hyperebenen so aufgeteilt, daß eine
Anzahl istr von Untermengen entsteht, so daß gilt: C(l) = {(i,j,r) ∊ G
| mod(i + j + k – 3,
istr) + 1 = l} (17)
-
Anders gesagt, die Untermenge C(l)
ist eine Vereinigung der Hyperebenen H(l), H(l+istr), H(l+2istr)...
-
Dann wird in der Reihenfolge C(1),
C(2),... C(istr), C(1), C(2),...C(istr) für die Vorwärtssubstitution und in der
Reihenfolge C(istr), C(istr-1),...C(1), C(istr), C(istr-1)...C(1) für die Rückwärtssubstitution jede
der Untermengen, die auch Farben (colors) genannt werden, berechnet.
Dies ist beispielhaft in 4 dargestellt.
Unbekannte Startwerte werden jeweils auf Null gesetzt. Bei diesem
Verfahren wird jede Farb-Untermenge zumindest zweimal berechnet,
so daß die
Abhängigkeit
von H(max(l-istr, 1)) von H(max(l-1), 1)) berücksichtigt wird.
-
Ein Beispiel für einen Algorithmus (2), der dieses
Verfahren durchführt,
ist:
-
(Vorwärtssubstitution)
-
- pi := 0 für
alle i
- für
id = 1, ..., nd; mit nd ≥ 2
- für
l= 1, ..., istr,
- Für
k = 1, ..., nzr
- Für
i = l, l + istr, l + 2 istr, ..., m
- ti := 0
- Ende der Schleife für
i
- Für
l, l + istr, l + 2 istrFür
i = , ..., m
- ti := ti + Li,i-izr_k pi-izr_k
- Ende der Schleife für
i
- Ende der Schleife für
k
- l + 2 istrFür
i = l, l + istr, , ..., m
- pi = (1/Δi,i)
(ri + ti)
- Ende der Schleife für
i
- Ende der Schleife für
l
- Ende der Schleife für
id
- und
-
(Rückwärtssubstitution)
-
- qi := 0 für alle i
- für
id = 1, ..., nd; mit nd >=
2
- für
l = istr, ..., 1,
- für
k = l, ..., nzr
- Für
i = l, l + istr, l + 2 istr, ..., m
- ti := 0
- Ende der Schleife für
i
- Für
i = l, l + istr, l + 2 istr, ..., m
- ti := ti + Ui,i+izr_k qi+izr_k
- Ende der Schleife für
i
- Ende der Schleife für
k
- Für
i = l, l + istr, l + 2 istr, ..., m
- qi = pi + (1/Δi,i)
ti
- Ende der Schleife für
i
- Ende der Schleife für
l
- Ende der Schleife für
id.
-
Ein Flußdiagramm entsprechend diesem Agorithmus
ist in 5 gezeigt.
-
Durch geeignete Wahl von istr kann
erreicht werden, daß jede
Farb-Untermenge soviele Gitterpunkte enthält, daß eine vektorielle Verarbeitung sinnvoll
möglich
ist. Für
den bekannten 7-Punkt-Stencil können
die Gitterpunkte darüberhinaus
so angeordnet weiden, daß für jedes
Element aus C(l) mit einem konstanten Abstand (stride) auf den Speicher
zugegriffen und somit die indirekte Adressierung vermieden werden
kann. Beispielsweise in dem Fall nx = ny und beliebigem nz kann
istr = nx-1 gewählt werden. In diesem Fall
können
alle Speicherzugriffe zur Berechnung der Gitterpunkte innerhalb
C(l) mit konstantem Abstand istr durchgeführt werden, und andererseits
ist (nx ny nz) / nc als Vektorlänge hinreichend groß.
-
Zum Beispiel kann im Fall des 7-Punkt-Stencils
mit einem Gitter, das in x-Richtung nx=10,
in y-Richtung ny=11 und in z-Richtung nz=10 Gitterpunkte hat, und mit dem wie oben
bestimmten istr = 9 eine Aufteilung in Farb-Untermengen vorgenommen
werden, die dann, wie in 4 gezeigt,
berechnet werden. Es ist somit klar, daß die zur Berechnung eines Gitterpunktes
i benötigten
6 weiteren Gitterpunkte, nämlich
i+1, i+10, i+110, i-1, i-10 und i-110 nicht derselben Farb-Untermenge angehören, und
daher ist die Berechnung eines Gitterpunktes einer Farb-Untermenge
unabhängig
von den Berechnungen der anderen Gitterpunkte der selben Farb-Untermenge. Eine
indirekte Adressierung ist ebenfalls nicht nötig, da allein aus der Kenntnis
von i sofort ersichtlich ist, welche weiteren Gitterpunkte benötigt werden.
Eine schnelle Bearbeitung auf einem Vektorprozessor ist möglich, da
die entsprechenden Speicherbereiche, in denen die Daten, der Gitterpunkte
i+1, i+10, i+110, i-1, i-10 und i-110 gespeichert sind, bereits
vorgeladen werden können.
-
Eine weitere Verbesserung dieser
Multi-Color-Technik, die in dem obengenannten Artikel ebenfalls
offenbart ist, ist die sog. "Overlapping
Multicolor"-Technik.
Hier wird die Konvergenz der Multi-Color-Technik dadurch verbessert,
daß für die Vorwärts-Substitution
die Farb-Untermengen in der Reihenfolge C(istr-w), (C(istr-w-1)...C(istr),
C(1), C(2) ... C(istr) berechnet werden. Die Überlappung, d.h. die Zahl der
Farb-Untermengen, die zusätzlich
berechnet werden, beträgt
w. Entsprechendes gilt für
die Rückwärtssubstitution.
Diese überlappende
Multi-Color-Technik zeigt eine verbesserte Konvergenz gegenüber der
einfachen Multi-Color-Technik,
bietet aber die gleichen Vorteile bezüglich Vektorlänge und Speicherzugriff.
-
Eine Schwierigkeit beider Multi-Color-Technik
und der Multi-Color-Technik mit Überlappung
ist die Bestimmung der Anzahl von Farb-Untermengen und die Anwendung
auf andere Fälle
als die bekannten 7-Punkt-Stencil oder 5-Punkt-Stencil. Eine zu kleine
Anzahl Farb-Untermengen führt
zu einem erhöhten
Rechenaufwand, da dann der Startwert für die erste Näherung zu
weit von. der Lösung
entfernt ist. Eine zu große
Anzahl Farb-Untermengen führt
zu einer geringeren Rechengeschwindigkeit, da die Vektorlängen für die Warteschlangenstruktur
des Vektorrechners zu kurz werden.
-
Die Bestimmung der Größe istr
ist bei anderen als den genannten 5-Punkt- und 7-Punkt-Stencil nicht
so leicht möglich,
da Voraussetzung bleibt, daß die
einzelnen Gitterpunkte einer Farbuntermenge unabhängig voneinander
sein müssen.
-
Aufgabe der Erfindung ist es, ein
verbessertes Verfahren zur Lösung
der Gleichung A q
= r (1)
auf
einem Vektorcomputer mit gegebener Warteschlangen-Länge zu schaffen.
-
Diese Aufgabe wird durch ein Verfahren nach
Anspruch 1 gelöst.
Die abhängigen
Ansprüche betreffen
weitere vorteilhafte Aspekte der Erfindung.
-
Das erfindungsgemäße Verfahren ist insbesondere
bei der Lösung
partieller elliptischer Differentialgleichungen für ein regelmäßiges quadratisches
Gitter geeignet.
-
Ein Vorteil der Erfindung ist, daß das Verfahren
sich speziell an die Gegebenheiten eines Vektorrechners, wie etwa
die Warteschlangen-Länge
und die Zahl der Speicher-Bänke,
anpaßt,
so daß es schneller
arbeitet als die bekannten Verfahren.
-
Ein weiterer Vorteil des erfindungsgemäßen Verfahrens
ist, daß ohne
weiteren Aufwand auch andere Abhängigkeiten
als die bekannten 5-Punkt-Stencil oder 7-Punkt-Stencil berücksichtigt werden
kann.
-
Noch ein weiterer Vorteil des erfindungsgemäßen Verfahrens
ist es, daß es
gut parallelisiert werden kann, d.h. auf mehreren Prozessoren parallel verarbeitet
werden kann, insbesondere auch in einem Computersystem mit verteiltem
Speicher.
-
Die Aufgabe, Vorteile und besondere
Aspekte der Erfindung werden anhand der Beschreibung einer speziellen
Ausführungsform
im Anschluß erklärt. In den
Zeichnungen zeigt:
-
1 eine
Ansicht eines zweidimensionalen Gitters mit einem 5-Punkt-Stencils;
-
2 eine
Ansicht eines dreidimensionalen Gitters mit einem 7-Punkt-Stencil,
-
3A eine
schematische Ansicht einer Matrix M für den Fall des 7-Punkt-Stencils;
-
3B eine
schematische Ansicht der Lage der Hyperebenen im dreidimensionalen
Raum;
-
4 eine
schematische Ansicht der Reihenfolge der Berechnung für den Fall
nx=10, ny=11 und
istr=9 bei beliebigem nz;
-
5 ein
Flußdiagramm
zur weiteren Erklärung
von Algorithmus 2;
-
6 ein
Basisdiagramm zur Erklärung
des Betriebs eines Vektorrechners;
-
7 eine
Kurve, die den Zusammenhang zwischen der Rechenleistung und der
Vektorlänge bei
einem bekannten Vektorrechner darstellt;
-
8 ein
Flußdiagramm
zur weiteren Erklärung
von Algorithmus 4;
-
9 eine
schematische Ansicht der Reihenfolge der Berechnungen für den Fall
nx=10, ny=11 und
istr=9 bei beliebigem nz und Verteilung
auf zwei parallel arbeitende Prozessoren;
-
10 ein
Flußdiagramm
zur weiteren Erklärung
von Algorithmus 5; und
-
11 und 12 weitere im dreidimensionalen Raum
mögliche
Stencils;
-
Beschreibung
der bevorzugten Ausführungsformen
-
Um das Verständnis der Erfindung zu erleichtern,
werden zunächst
einige Grundlagen von Vektorrechnern erläutert.
-
Bekannte Computer enthalten einen
oder mehrere Prozessoren, einen oder mehrere Speicher, Eingabe-/Ausgabevorrichtungen
und Kommunikationskanäle
zwischen diesen Bauteilen.
-
Um einen Befehl auszuführen, muß der Prozessor
diesen Befehl aus einem Speicher laden, dekodieren und ausführen. Die
für einen
solchen Zyklus benötigte
Zeit hängt – neben
der Komplexität
des Befehls und der internen Bauweise des Prozessors – auch von
der Organisation der Speicher und der Kommunikationskanäle ab.
-
Computerprogramme durchlaufen oft
Schleifen, die nur einen oder wenige einfache Befehle enthalten.
Ein Beispiel einer solchen Schleife ist die Berechnung einer Vektorsumme
und/oder eines Vektorkreuzproduktes.
-
Hierbei ist zu beachten, daß auch die
einfache Multiplikation zweier Zahlen in einem binär arbeitenden
Prozessor für
gewöhnlich
eine Reihe einzelner Schritte umfaßt, z.B. im einfachsten Fall
das Verschieben einer Bitfolge nach links um ein Bit (Multiplikation
mit 2) und die Addition binärer
Zahlen.
-
Wird ein solcher Algorithmus mit
einem einfachen. bekannten skalaren Computer bearbeitet, so wird
ein Großteil
der Zeit dafür
verbraucht, immer wieder den gleichen Befehl zu laden, zu dekodieren und
Zwischenergebnisse zu speichern.
-
Um dieses Problem zu lösen, werden
Vektorcomputer eingesetzt.
-
Die Grundstruktur des Betriebs eines
Vektorcomputers ist in 6 dargestellt.
-
In einem ersten Takt wird in dem
gezeigten Beispiel, ein erster Befehl aus einem Befehlsspeicher entnommen.
Im zweiten Takt wird dieser Befehl an jenen Teil des Prozessors
weitergereicht, der die Operanden bereitstellt, und ein zweiter
Befehl wird aus dem Befehlsspeicher entnommen. Vom dritten bis m-ten
Takt wird der erste Befehl ausgeführt, vom vierten bis (m + 1)-ten
Takt der zweite Befehl usw. Nach m-Takten gibt der Prozessor mit
jedem Takt ein Ergebnis aus, so daß für die Verarbeitung von n-Befehlen
insgesamt n + m – 1
Takte benötigt
werden. Ein skalarer Computer hingegen würde (n m) Takte benötigen.
-
Mit anderen Worten verwenden Vektorcomputer
Warteschlangen, wobei Daten/Befehle von einem Speicher in ein Vektorregister über eine
Lade-Warteschlange übertragen
werden. Die Berechnung wird dann in einer arithmetischen Warteschlange
ausgeführt,
und das Ergebnis wird zurück
in den Speicher über
Speicher-Warteschlangen geliefert. Die arithmetische Warteschlange
ist eine multifunktionelle Warteschlange, die einzelne Schritte – wie eine
Gleitkomma-Berechnung – parallel
ausführen kann.
-
Die Effizienz oder Rechenleistung
eines solchen Vektorcomputers ausgedrückt in Gleitkomma-Rechenoperationen
pro Sekunde (FLOPS) hängt stark
von der sog. Vektorlänge
ab. Aus 6 ist zu erkennen,
daß bei
dem vektoriell zu verarbeitenden Befehlsstapel der erste Befehl
m Takte benötigt,
während
jeder weitere Befehl nur noch einen zusätzlichen Takt hinzufügt. Diese
Verarbeitungweise setzt voraus, daß der n-te Befehl ausgeführt werden
kann, ohne daß das
Ergebnis eines der vorangehenden m-1 Befehle bekannt sein muß, wie das
typischerweise bei den obengenannten Vektorverknüpfungen der Fall ist.
-
Unter Vektorlänge wird deshalb eine für die Anzahl
vektoriell ausführbarer
Befehle repräsentative Zahl
verstanden. Üblicherweise
kann die geeignete Vektorlänge
für einen
gegebenen Vektorcomputer durch eine einfache Schleife abgeschätzt werden. Als
Beispiel kann folgender Algorithmus (3) dienen:
n := 2
MFlopsmax
:= 0
Beginn der Schleife über
n
Zurücksetzen
der inneren Uhr
Für
i = 1 bis n
xi=xi +
pi zi
Ende
der Schleife für
i
MFlops = n / Zeit
Wenn MFlops > MFlopsmax
MFlopsmax := MFlops
Sonst
lvc
:= n
Ende der Abfrage
n = n + 1
Ende der Schleife
für n.
-
7 zeigt
am Beispiel eines NEC FX4-Computers den Einfluß der Vektorlänge auf
die Rechenleistung. Beobachtet man zunächst einen deutlichen Anstieg
der Rechenleistung mit zunehmender Vektorlänge, so zeigt oberhalb von
etwa n = 400 die Vektorlänge
kaum noch Einfluß auf
die Rechenleistung.
-
In der vorliegenden Beschreibung
wird unter Mindest-Vektorlänge
jene Vektorlänge
verstanden, ab der die Rechenleistung annähernd unabhängig von der Vektorlänge ist.
In dem gezeigten Beispiel ist dies etwa 400.
-
Diese Mindest-Vektorlänge hängt von
dem jeweils verwendeten Vektorprozessor ab. So benötigt etwa
der Fujitsu-VPX Vektoren mit einer Mindestlänge von 1000, um die Vorteile
der Vektorisierung voll ausnutzen zu können, während der Cray-EL98 mit deutlich
kürzeren
Vektorlängen
auskommt.
-
Daß oberhalb der jeweiligen Mindest-Vektorlänge kaum
ein Einfluß der
Vektorlänge
auf die Rechenleistung zu beobach ten ist, ist auf eine Reihe von
Gründen
zurückzuführen. So
ist offensichtlich, daß bei
voller Überlappung
in 6, d.h. wenn die Zahl
der Stufen der Warteschlange gleich der Zahl der Befehle eines Stapels
ist, mit einer weiteren Zunahme der Vektorlänge nur noch ein geringer Vorteil erreicht
werden kann.
-
Eine weitere die Rechenleistung eines
Computers beeinträchtigende
Größe ist die
Anzahl, der Aufbau und die Organisation des Speichers. Um auf einen
Speicherplatz zuzugreifen, wird eine zugehörige Speicherbank für gewöhnlich vorgeladen,
dann erst erfolgt der Zugriff. Anschließend wird die Speicherbank
wieder heruntergefahren. Die Zeit zum Vorladen und Herunterfahren
kann ein Vielfaches der eigentlichen Zeit des Zugriffs betragen.
Deshalb ist es üblich,
daß die
Bits, auf die nacheinander zugegriffen werden soll, in verschiedenen
Speicherbänken
gespeichert werden, so daß immer
im Voraus die nächsten
Speicherbänke
vorgeladen werden können.
Hierdurch wird eine beachtliche Verkürzung der insgesamt benötigten Speicherzugriffszeit
erreicht. Allerdings. setzt dies voraus, daß immer bekannt ist, auf welche
Speicherbänke
demnächst
zugegriffen werden soll. Bei dem eingangs geschilderten Hyper-Plane-Verfahren
mit der indirekten Adressierung ist dies unmöglich, da die Adressen der
nächsten
Zugriffe immer erst berechnet werden müssen. Es gilt daher, indirekte
Adressierungen zu vermeiden.
-
Optimal hingegen ist ein Speicherzugriff
mit bekanntem und konstantem Stride, d.h. also mit einem konstanten,
immer schon bekannten Abstand zwischen nacheinander aufgerufenen
Speicherplätzen.
-
Ähnliche
Speicherzugriffsprobleme treten auf, wenn trotz direkter Adressierung
kurz hintereinander auf diesselbe Speicherbank zugegriffen werden
soll. Auch dies vermindert die Rechenleistung, da die Speicherbank
vor einem erneuten Zugriff erst heruntergefahren und dann wieder
vorgeladen werden muß.
-
Vektoriell zu verarbeitende Schleifen
werden deshalb vorzugsweise so ausgebildet, daß zwischen zwei Zugriffen auf
eine Speicherbank hinreichend Zeit vergeht.
-
Eine weitere Verbesserung der Rechenleistung
kann erzielt werden, wenn die Bearbeitung der Befehle auf mehrere
unabhängige
parallel arbeitende Vektorprozessoren verteilt wird. Hierfür ist Voraussetzung,
daß mehrere
unabhängige
Teilaufgaben gebildet werden können.
-
Wie eingangs geschildert wurde, betrifft
die Erfindung die Lösung
der Gleichung 1 mit Hilfe der Multi-Color-Technik auf einem Vektorcomputer.
-
Um den vorangehend aufgezeigten Einflußgrößen auf
die Rechenleistung gerecht zu werden, schlägt die Erfindung vor, einen
konstanten Speicher-Zugriffsabstand istr nach folgenden Regeln zu bestimmen:
- a) mod(izrk, istr)≠0 für alle k=1...nzr;
mit izrk= {j:Li,i-j≠0 oder Ui,i+j≠0;}
und wobei nzr die Anzahl aller j ist für die gilt, Li,i-j≠0 oder Ui,i+j≠0;
- b) der größte gemeinsame
Teiler von istr und nb ist 1 wobei nb die Anzahl der Speicherbänke ist;
- c) istr ist möglichst
klein, so daß die
tatsächliche Vektorlänge größer als
die für
den Vektorprzessor benötigte
Mindestvektorlänge
ist; und
- d) der tatsächliche
Wert für
istr ergibt ein Maximum für
die Größe x, wobei
gilt: x = min{int(istr/mod(izrk, istr))
mit k=1,...nzr}
-
Die Regel a) stellt sicher, daß alle Rechenoperationen
innerhalb einer Farb-Untermenge unabhängig voneinander ausge führt werden
können,
so daß überhaupt
innerhalb einer Farbuntermenge eine vektorielle Verarbeitung möglich wird.
Dies setzt voraus, daß keine
der Zahlen izrk ohne Rest durch istr teilbar
ist. Durch Anwenden dieser Regel ist es möglich, nicht nur den 5-Punkt-Stencil
oder den 7-Punkt-Stencil zu handhaben, sondern beliebige Abhängigkeiten
können
durch das erfindungsgemäße Verfahren
berücksichtigt
werden.
-
Die Regel b) stellt sicher, daß kein Speicherkonflikt
auftritt, also daß nicht
hintereinander auf diesselbe Speicherbank zugegriffen wird und daß alle Speicherbänke ausgenutzt
werden.
-
Die Regel c) stellt sicher, daß die Vektorlänge groß genug
ist, um den Vektorcomputer optimal zu nutzen. Hierbei ist jedoch
zu beachten, daß mit
steigender Vektorlänge
die Konvergenz des Verfahrens verschlechtert wird. Dies ist darauf
zurückzuführen, daß zu Beginn
sowohl der Vorwärts-
als auch der Rückwärts-Substitution
die Größen pi bzw. qi für alle i jeweils
auf den Wert 0 (oder einen anderen geeigneten Startwert) gesetzt
sind. Nach dem ersten Durchlaufen der Schleife über l sind dann die Größen pl, pl+istr, pl+2istr,... pm bzw.
ql, ql+istr, ql+2istr,....qm jeweils auf einem neuen Wert,
dem jeweiligen entsprechenden ersten Näherungswert, der dann in alle
weiteren Berechnungen der übrigen
Größen pi bzw. qi eingeht.
Die Konvergenz ist umso effektiver, je schneller die ersten Näherungen
für die
weitere Berechnung zur Verfügung
stehen.
-
Um eine gute Konvergenz sicherzustellen, ist
deshalb die Regel d) vorgesehen. Diese Regel ist mit der geschilderten
Von-Neumann-Reihenentwicklung verknüpft und stellt sicher, daß die Konvergenz des
erfindungsgemäßen Verfahrens
zumindest so gut wie eine Reihenentwicklung nach Gleichung (16) bis
zum Glied ntr = (nd – 1)
x ist.
-
Ein bevorzugter Algorithmus (4) zur
Berechnung von istr ist:
istrmax := int(m / lvc)
x :=
0
istr := 0
für
i = 2, 3, ..., istrmax
Wenn ggT(i, nb) = 1, dann
modmin
:= i
modmax := 0
für
k = 1, ..., nzr
modmin := min modmin, mod(izrk,
i))
modmax := max(modmax, mod(izrk,
i))
Ende der Schleife für
k
Wenn modmin ≢ 0,
dann
xtest = int(i / modmax)
Wenn xtest > x, dann
x :=
xtest
istr := i
Ende der Abfrage xtest
Ende der Abfrage
ggT(i, nb)
Ende der Abfrage modmin
Ende der Schleife für i.
-
Ein Flußdiagramm entsprechend diesem
Algorithmus ist in 8 zur
weiteren Erklärung
gezeigt.
-
Hierbei ist lvc die Mindest-Vektorlänge, also die
minimale für
den Vektorprozessor zu fordernde Vektorlänge, und ggT(i, nb) der größte gemeinsame Teiler
von i und nb.
-
Alternativ dazu können insbesondere die letzten
beiden Regeln über
eine Fuzzy-Verknüpfung berücksichtigt
werden, da. die Mindest-Vektorlänge keine
scharfe Grenze ist. Im obengenannten Algorithmus (4) kann dies berücksichtigt
werden, indem beispielsweise die Schleife über i nicht bis istrmax, sondern
bis zu einem beliebigen Vielfachen von istrmax läuft, und wobei Lösungen mit
i > istrmax umso weniger
bevorzugt werden, als i größer als istrmax
ist.
-
Im Anschluß wird die Anwendung der oben beschriebenen
Erfindung im Zusammenhang mit einem Parallel-Vektorrechner mit np
= 2 Vektorprozessoren beschrieben. Für den Fachmann ist es offensichtlich,
daß in
analoger Art auch mehr als zwei Prozessoren ausgenützt werden
können.
-
Ein Beispiel mit zwei Prozessoren
und wiederum dem 7-Punkt-Stencil
beim Gitter aus 2 mit nx
= 10, ny = 11 und istr = 9 bei beliebigem nz ist in 9 gezeigt.
-
Zunächst berechnet der erste Prozessor
unabhängig
vom zweiten Prozessor die ersten Näherungslösungen für die Farb-Untermengen C(1)
und C(2). Gleichzeitig berechnet der zweite Prozessor unabhängig vom
ersten Prozessor die ersten Näherungslösungen für die Farb-Untermengen
C(5) und C(6).
-
Hierzu werden die an sich benötigten Werte für die Gitterpunkte
der Untermengen C(8), C(9) bzw. C(3), C(4) jeweils auf Null gesetzt.
-
Nach der Berechnung der Untermengen C(2)
bzw. C(6) findet eine Synchronisation statt. Das heißt, mit
der Berechnung der nächsten
Untermengen C(3) bzw: C(7) wird erst begonnen, wenn beide Prozessoren
ihre jeweiligen Berechnungen der vorangehenden Untermengen beendet
haben. Dies ist nötig,
da die Werte der Untermengen C(3) bzw. C(7) in die Berechnung der
Werte für
C (5) bzw . C (1) eingehen . Eine Änderung der Werte für C(3),
während die
Berechnung für
C(5) noch nicht abgeschlossen ist, würde zu einer Unsicherheit in
der Berechnung führen,
da die Werte aus C(3) für
die Berechnung von C(5) benötigt
werden, so daß das
Ergebnis für
C(5) dann von mehr oder minder zufälligen Arbeitsbedingungen der
jeweiligen Prozessoren abhängig
wäre.
-
Nach der Synchronisation werden dann
die Untermengen C(3) und C(4) durch den ersten Prozessor und die
Untermengen C(7), C(8) und C(9) durch den anderen Prozessor berechnet.
Anschließend
werden die Prozessoren erneut synchronisiert, und das Verfahren
beginnt mit einem neuen Durchlauf und der Berechnung der Untermengen
C(1) bzw. C(5), wobei jetzt die ersten Näherungslösungen für die Untermengen C(8), C(9)
bzw. C(3), C(4), die im vorangehenden Durchlauf berechnet wurden,
eingesetzt werden. Diese Schleife wird solange wiederholt, bis die
Bedingung id = nd in Algorithmus (2) erreicht ist, also bis davon
ausgegangen werden kann, daß die
Näherungslösung hinreichend.
gut ist. Üblicherweise
reichen zwei Durchläufe
(nd = 2) für
eine zufriedenstellende Näherungslösung aus.
-
Ein bevorzugter Algorithmus (5) für diese
Bearbeitung auf einem Paralell-Vektor-Rechner ist:
Für i= int((ip-1)
m / np) + 1, ..., int((ip – m)
/ np)
pi := 0
Ende der Schleife
für i
-
SYNCHRONISATION
-
(Vorwärtssubstitution)
-
- Für
id = 1, nd
- Für
ia = 0, 1
- Für
l = int((2 (ip-1) + ia) istr / (2 np)) + 1, ..., int((2 (ip-1) +
ia + 1) istr / (2 np))
- Für
i = l, l + istr, 2 + 2 istr, ..., m
- ti := 0
- Ende der Schleife für
i
- Für
k = 1, ..., nzr
- Für
i = l, l + istr, l + 2 istr,..., m
- ti := ti + Li,i-izr_k pi-izr_k
- Ende der Schleife für
i
- Ende der Schleife für
k
- Für
i = l, l + istr, l + 2 istr, ..., m
- pi = (1/Δi,i)
(ri + ti)
- Ende der Schleife für
i
- Ende der Schleife für
1
-
SYNCHRONISATION
-
- Ende der Schleife für
ia
- Ende der Schleife für
id
-
(Rückwärtssubstitution)
-
Für
i= int((ip-1) m / np) + 1, ..., int((ip – m) / np)
-
qi := 0
-
Ende der Schleife für i
-
SYNCHRONISATION
-
- Für
id = 1, ..., nd
- Für
ia = 1, 0
- Für
l = int((2 (ip-1) + ia) istr / (2 np)) + 1, ..., int((2 (ip-1) +
ia + 1) istr / (2 np)):
- Schritt -1
- ti := 0
- Ende der Schleife für
i
- Für
k = 1, ..., nzr
- Für
i = l, l + istr, l + 2 istr, ..., m
- ti := ti + Ui,i+izr_k qi+izr_k
- Ende der Schleife für
i
- Ende der Schleife für
k
- Für
i = 1, 1 + istr, 1 + 2 istr, ..., m
- qi = pi + (1/∆i,i)
ti
- Ende der Schleife für
i
- Ende der Schleife für
l
-
SYNCHRONISATION
-
- Ende der Schleife für
ia
- Ende der Schleife für
id.
-
Ein Flußdiagramm entsprechend diesem
Algorithmus, ist in 10 zur
weiteren Erklärung
gezeigt.
-
Wenn mehr als 2 Prozessoren vorgesehen sind,
laufen die Schleifen über
ia von 0 bis int((istr / np)) – 1,
wobei np die Anzahl der Prozessoren ist.
-
Der Fachmann erkennt leicht, daß das erfindungsgemäße Verfahren
sich an spezielle Gegebenheiten eines jeweiligen Vektorcomputers
anpaßt,
wie etwa die Warteschlangen-Länge
und die Zahl der Speicherbänke,
so daß es
schneller arbeitet als die bekannten Verfahren. Auch können ohne
weiteren Aufwand andere Abhängigkeiten
als die bekannten 5-Punkt-Stencil
bzw. 7-Punkt-Stencil berücksichtigt werden.
-
Beispielsweise können die in den 11 und 12 gezeigten 19- und 27-Punkt-Stencil
berechnet werden.
-
Desweiteren ist das erfindungsgemäße Verfahren
gut für
die parallele Verarbeitung auf mehreren unabhängigen Prozessoren geeignet,
insbesondere auch bei einem Computersystem mit verteiltem Speicher,
wobei jedem Prozessor ein eigener Speicherbereich zugewiesen ist,
so daß Zugriffskonflikte
vermieden werden.