lösung linearer gleichungssysteme · pdf file1 1 einleitung die lösung linearer...
Post on 06-Feb-2018
225 Views
Preview:
TRANSCRIPT
Westfälische Wilhelms-Universität Münster
Ausarbeitung
Lösung linearer Gleichungssysteme
(im Rahmen des Seminar „Parallele Programmierung“)
Nico Ziborius
Themensteller: Prof. Dr. Herbert KuchenBetreuer: (Prof. Dr. Herbert Kuchen)Institut für WirtschaftsinformatikPraktische Informatik in der Wirtschaft
Westfälische Wilhelms-Universität Münster
III
Inhaltsverzeichnis
1 Einleitung................................................................................................................... 1
2 Direkte Verfahren ...................................................................................................... 2
2.1 Gauß-Elimination.................................................................................................. 2
2.2 zyklische Reduktion.............................................................................................. 7
3 Iterative Verfahren................................................................................................... 13
3.1 Klassische iterative Verfahren ............................................................................ 13
3.1.1 Jacobi-Verfahren ............................................................................................ 143.1.2 Gauß-Seidel-Verfahren .................................................................................. 16
3.2 Methode der konjugierten Gradienten ................................................................ 20
4 Zusammenfassung ................................................................................................... 23
A Anhang Implementierung der Gauß-Elimination .................................................... 25
Literaturverzeichnis ........................................................................................................ 27
1
1 Einleitung
Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer
Anwendungen. Viele naturwissenschaftliche Phänomene werden durch partielle
Differentialgleichungen beschrieben werden, deren numerische Lösung oft auf lineare
Gleichungssysteme zurückzuführen ist. Große Gleichungssysteme treten z. B. in den
Ingenieurswissenschaften bei der Simulation von Flüssigkeits- oder Gasströmungen
oder bei der Simulation des Verhaltens von Materialien wie Stahl, Beton, Holz usw.
unter Belastung auf.
Gesucht wird jeweils die Lösung eines linearen Gleichungssystems Ax =b wobei A eine
nichtsinguläre n×n Matrix darstellt, b einen aus n Einträgen bestehenden Vektor und x
den aus n Einträgen bestehenden Lösungsvektor bezeichnet. Die Lösungsverfahren
lassen sich in zwei Klassen, Direkte-Verfahren und Iterative-Verfahren, einteilen.
Direkte Verfahren liefern, in einer von der Größe des Gleichungssystems abhängigen
Anzahl von Berechnungsschritten, eine exakte Lösung. Iterative-Verfahren berechnen
hingegen eine Näherung an die Lösung, in dem sie, ausgehend von einem frei zu
wählenden Startvektor, eine Folge von Vektoren berechnen die gegen die gesuchte
Lösung konvergiert.
In Kapitel eins werden zwei direkte Verfahren beschrieben, eine parallele
Implementierung des Gauß-Eliminations-Verfahrens und die zyklische Reduktion, die
sich für die Lösung linearer Gleichungssysteme eignet, deren Matrix A eine
Bandstruktur aufweist. In Kapitel 3. werden die klassischen iterativen Verfahren Jacobi-
Verfahren und Gauß-Seidel-Verfahren beschrieben, sowie die Methode der konjugierten
Gradienten. Die parallele Implementierung der vorgestellten Algorithmen ist jeweils für
Rechner mit verteiltem Adressraum ausgelegt.
2
2 Direkte Verfahren
2.1 Gauß-Elimination
Gesucht wird die Lösung für das Gleichungssystem Ax =b mit der nichtsingulären
Matrix A∈Rn×n und Vektor b∈Rn . Idee der Gauß-Elimination ist es durch
Transformation das Gleichungssystems Ax = b in ein äquivalentes Gleichungssystem
A‘x = b‘ zu überführen, so das A‘ eine obere Dreiecksmatrix bildet.
Die in Schritt k der Transformation entstehende Matrix zeichnet sich dadurch aus, dass
in den ersten k-1 Spalten unterhalb der Diagonalelemente nur Nullen stehen.
Matrix A(k+1) und Vektor b(k+1) werden aus A(k) , b(k) gebildet in dem geeignete
Vielfache der k-ten Zeile von A(k) bzw. des k-ten Elementes von b(k) zu den Zeilen
k+1, ..., n addiert werden, so das die Koeffizienten von xk in den Gleichungen k+1, ..., n
Null gesetzt werden. Hierzu werden zunächst Eliminationsfaktoren
ermittelt.
�������
�
�
�������
�
�
=
a
aaaaa
nn
n
n
A
000
00
´222
12111
�
���
�
nikfür )(
)(
≤<=aal k
kk
k
ikik
�����������
�
�
�����������
�
�
=−
−
−
−
−
−−
−
−
aa
aaaaa
aaaaaaaaa
k
nn
k
nk
k
kn
k
kk
k
nk
k
kk
k
kk
nkk
nkk
k
)()(
)()(
)1(
,1
)1(
,1
)1(
1,1
)2(
2
)2(
2
)2(
1,2
)2(
22
111,11211
)(
00
0
0
A
��
�����
��
��
������
�
�
Abbildung 1 obere Dreiecksmatrix
Abbildung 2 Matrix A nach
Transformationsschritt k
Gleichung 1
3
��
�
�
�−=
+=
n
kj
n
kj
n
kn
kk
k aba 1
j
)()(
)( xx 1
Danach werden die Elemente der Matrix A sowie des Vektors b neu berechnet
für k< j ≤ n und k <i ≤ n. All Zeilen oberhalb der k-ten Zeile bleiben somit unverändert.
Nach n-1 Schritten erhält man die obere Dreiecksmatrix A(n). Durch
Rückwärtseinsetzten können nun xn, xn-1, ..., x1 nacheinander bestimmt werden, hierzu
wird für k = n,n-1,...,1
ermittelt.
Die Gauß-Elimination funktioniert also nur, wenn die Elemente auf der Diagonalen
ungleich Null sind, da sonst Division durch Null in den Gleichungen 1 und 4 zu einem
Fehler führen würde, obwohl das Gleichungssystem prinzipiell lösbar ist. Sind die
Elemente auf der Diagonalen der Matrix A sehr klein, so werden die daraus
resultierenden Eliminationsfaktoren dementsprechend sehr groß, was zu Rundungs-
fehlern und somit zu Ungenauigkeiten, führen kann. Daher sind Strategien erforderlich
um die Elemente auf der Diagonalen der Matrix durch ein andere, sogenanntes
Pivotelement, zu ersetzten. Hierfür sind verschiedene Strategien denkbar, die
Spaltenpivotsuche sucht aus den Elementen akk(k),..,ank
(k) das betragsmäßig größte ark(k)
aus und vertauscht dann die Zeilen k und r , falls r ≠ k. Andere Strategien sind die
Zeilenpivotsuche, die analog zur Spaltenpivotsuche arbeitet, oder die Totalstrategie,
welche das betragsmäßig größte Element innerhalb der Teilmatrix Ã(k) = (aij(k)) mit k≤
i,j ≤ n sucht.
Eine sequentielle Implementierung des Gauß-Eliminations-Verfahrens benötigt zur
Berechnung der oberen Dreiecksmatrix drei ineinander verschachtelte Schleifen, deren
innere Schleife ungefähr n3/3 mal durchlaufen wird. Laufzeit dieser Implementierung
läßt sich somit zu O(n3) abschätzten, eine sequentielle Implementierung ist z.B. in
[RA98] zu finden.
alaak
kjik
k
ij
k
ij
)()()1(*−=
blbbk
kik
k
i
k
i
)()()1(*−=
+
Gleichung 2
Gleichung 3
Gleichung 4
4
Eine Möglichkeit der parallelen Implementierung diesem Algorithmus erfolgt durch,
zeilenblockweise Aufteilung der Matrix A auf die einzelnen Prozessoren. Sei p die
Anzahl der Prozessoren und n durch p ganzzahlig Teilbar, dann erhält Prozessor Pi für 1
≤ i ≤ p, eine n/p großen Zeilenblock zugewiesen. Jeder Prozessor Pi hat ein 2-
dimensionales Array a[n/p][n] mit den Koeffizienten seines Teils der Matrix A, ein
Array b[n/p] mit den entsprechenden Werten der Vektors b, ein Array marked[n/p], um
sich zu merken ob Zeile r schon Pivot Zeile war und ein weiters Array pivot[n/p], in
dem festgehalten wird in welchem Transformationsschritt eine Zeile Pivotzeile war. Zur
Lösung des Gleichungssystems wird zunächst das Array marked mit 0 initialisiert.
Danach wird in n-1 Schritten die obere Dreiecksmatrix berechnet. Zur Aufstellung der
Zwischenmatrix A(k+1) aus A(k) sind dabei jeweils folgende Berechnungen, von jedem
Prozessor Pi, vorzunehmen:
1. lokales Pivotelement bestimmen: Pi ermittelt unter den Elementen, akk(k),..,ank
(k) ,
die er von Spalte k hält, lokal das Maximum, dabei werden die Werte ark(k) nicht
berücksichtigt, wenn Zeile r schon Pivotzeile war, also marked[r]=1.
2. globales Pivotelement bestimmen: Aus den lokalen Maxima aus Schritt 1. wird
das globale Maximum ark(k) ermittelt, indem Pi die Funktion
MAX_TOURNAMENT(int i, int value, int winner) aufruft. Die Variablen i, winner
enthalten dabei die eigene Prozessor-ID und der Parameter value das lokale
Maximum. MAX_TOURNAMENT in Pseudocode:
int MAX_TOURNAMENT(int i, int value,int winner){int j,k;for(k=0; k<log p-1; k++){
j=i^(1<<k);[j]tmp_value ← value;[j]tmp_winner ← winner;if(tmp_value>value){ value = tmp_value; winner = tmp_winner;}
}return winner;
}
1
0
0
2
4
-1 -5 11
4
3 -4
-10
-8 -2 10
7
-6
14
3
0
1
0
0
-
1
-
-
ba
Abbildung 3 Beispiel mit 4×4 Matrix A
aufgeteilt auf 2 Prozessoren in Schritt k=2 hat P1
als lokales Maximum 4 aus Zeile 2, Zeile 1 wird
nicht berücksichtigt, da sie Pivotzeile in Schritt
k=1 war. P2 ermittelt aus |–8| und |–1|, |-8| als
lokales Maximum
5
MAX_TOURNAMENT implementiert den „Turnier“- Algorithmus, Rückgabewert
ist die ID des Prozessors mit dem größten Wert. In jedem Schritt enthält value den
aktuell größten Wert und winner die dazugehörige Prozessor-ID. Austausch der
Werte findet in einem „two message-passing step“ statt.
3. Markieren und Verteilen der Pivotzeile: Falls Pi der Gewinner aus 2. und ark(k)
sein lokales Maximum, dann markiert er Zeile r als Pivotzeile, in dem er
marked[r]=1 und pivot[r]=k setzt, kopiert er die Werte ark, ..., arn und b[r] in einen
Hilfsvektor und sendet diesen mittels einer Broadcastoperation an alle andere
Prozessoren, damit diese die nachfolgenden Berechnungen durchführen können.
4. Berechnung der Eliminationsfaktoren: Pi berechnet nun für alle seine lokalen
Zeilen i, mit i =1,...n/p, die noch nicht Pivotzeile waren (also marked[i]=0), lokal
die Eliminationsfaktoren lik nach Gleichung 1.
5. Neu Berechnung Matrixelemente: Pi berechnet für seine lokal gehaltenen
Elemente von A(k) und b(k), mit Hilfe der Eliminationsfaktoren, seine Werte der
Matrix A(k+1) und des Vektors b(k+1) laut den Gleichungen 2 und 3. Hierbei werden
nur die Werte aus den Zeilen neu berechnet, die noch nicht Pivotzeile waren.
Im Array pivot ist nun für alle Zeilen eingetragen, wann sie Pivotzeile waren, nur eine
Zeile war nie Pivotzeile, daher sucht jeder Prozessor in seinem Array marked nach der
1
0
0
2
4
-1 -5 11
4
3 -4
-10
-8 -2 10
7
-6
14
3
0
1
0
1
-
1
-
2
ba
1
0
0
2
4
-1 -5 11
4
3 -4
-10
-8 -2 10
7
-6
14
3
0
1
0
1
-
1
-
2
ba
Abbildung 4 Fortsetzung des Beispiels aus
Abb. 3. Nach dem 8 aus Zeile 3 als globales
Maximum ermittelt wurde setzt P2 für diese
Zeile marked=1 und pivot=2
Abbildung 5 Fortsetzung des Beispiels aus
Abb. 3. P1 und P2 müssen jeweils einen
Eliminationsfaktor für Zeile 2 bzw. 4 berechnen
P1 rechnet (4/-8), P2(-1/-8). In 5. muß
abschließend P1, P2 jeweils 2 Werte, in Zeile 2
und 4 neu berechnen.
6
unmarkierten Zeile j und setzt für diese pivot[j]=n. Eine Implementierung in
Pseudocode ist in Anhang A aufgeführt.
Nachdem so Matrix A in eine obere Dreiecksmatrix überführt wurde, kann die Lösung
des Gleichungssystems durch Rückwärtseinsetzten ermittelt werden, wofür wiederum n
Schritte notwendig sind. In jedem Schritt j berechnet der Prozessor der Zeile j enthält xj
gemäß Gleichung 4 und sendet sein Ergebnis an alle anderen Prozessoren die
Reihenfolge der Zeilen j wird dabei durch das Array pivot festgelegt, begonnen wird
mit der Zeile j für die gilt pivot[j]=n, dann folgt n-1,..,1. Aufgrund der
Datenabhängigkeiten, die beim Rückwärtseinsetzten gegeben sind, kommt es in dieser
Phase zu einer Sequentialisierung des Algorithmus. Eine Implementierung in
Pseudocode ist in Anhang A aufgeführt.
Als Rechenaufwand für die Bildung der oberen Dreiecksmatrix wird für die
Initialisierung des Array marked[n/p] ein Aufwand von � (n/p) benötigt, jeder der n-1
Iterationsschritte haten eine wird ein Aufwand für die Ermittlung des lokalen Maximuns
von � (n/p), � (n-1) für das Kopieren der Pivozeile und � [(n/p) ⋅(n-1)] für die
Berechung der Eliminationsfaktoren und die Neuberechnung der Matrixelemente, jeder
Prozessor hat n/p Zeilen und muss in jeder maximal n-1 Werte neu ermitteln. Insgesamt
stellt sich der Rechenaufwand für die Ermittlung der oberen Dreiecksmatrix wie folgt
dar:
Dies läßt sich umformen zu:
Damit ergibt sich eine fast Kostenoptimale Aufteilung auf die Prozessoren. Hinzu
kommt noch ein Kommunikationsaufwand von
.)1()(),(1
11 �
�
�
�
�+���
�
�−+−++Θ=
−
=
n
i pn
npn
inpn
pn
pnt
( ) ( )
��
�
�
�+Θ=
���
�
�+
−+
++Θ=
��
�
�
�+++Θ=
��
�
�
�+−−++Θ=
=
−
=
−
=
−
=
23
2
1
1
1
1
1
11
2)1(
2)1(
1),(
np
n
pnnnnn
pn
pn
pn
iipn
pn
pn
ininpn
pn
pntn
i
n
i
n
i
n
i
( ) ,loglog)(log),(1
12 �
�
�
�
�+−+Θ=
−
=
n
i
ppinppnt
7
der aus der Ermittlung des globalen Maximums zur Bestimmung der Pivotzeile und
dem Verteilen der Pivotzeile resultiert. Durch Umformung ergibt sich
Wegen t1<<t2 wird für große n der Rechenaufwand den Kommunikationsaufwand
übersteigen. Daher ist ein fast linearer Speedup möglich.
Vorteil dieser Implementierung der Gauß-Elimination gegenüber anderen, die eine
zeilenzyklische, oder gesamtzyklische Verteilung der Matrix A auf die Prozessoren
vornehmen ist, dass die Pivotzeile nicht zwischen zwei Prozessoren getauscht wird
sondern nur im Array marked, als solche gekennzeichnet. Eine zeilenzyklische bzw.
gesamtzyklische Implementierung ist in [RA98] zu finden. Durch diese
Vorgehensweise wird zusätzlicher Kommunikationsaufwand vermieden . Allgemein
liegen die Vorteile von parallelen Implementierungen der Gauß-Elimination in der
Vorhersagbarkeit der Laufzeit und des Speicherbedarfs, sowie der Berechnung einer,
abgesehen von Rundungsfehlern, exakten Lösung. Nacheilt ist, der insbesondere bei
dünnbesetzten Matrizen entstehende fill-in, durch Umformung der Matrix kann es zu
einer Auffüllung mit Nichtnullelementen kommen und so zu einer dichter besetzten
Matrix. Diese hat sowohl Auswirkungen auf den benötigten Speicherplatz, wie auch auf
den Rechenaufwand. Daher ist das Gauß-Elimination für dünnbesetzte Matrizen wenig
geeignet.
2.2 zyklische Reduktion
Sehr große Lineare Gleichungssystem weisen oftmals eine Bandstruktur auf, die
dadurch gekennzeichnet ist, dass alle Elementen, außer auf der Diagonalen und einigen
Nebendiagonalen, gleich Null sind. Eine Matrix A mit A= (aji) i,j=1,...,1 n ∈Rn×n heißt
Bandmatrix mit halber Bandbreite r , falls aij = 0 für |i-j| >r.
( ) ( )
( )pnpnn
pnppinppntn
i
n
i
loglog2
)1(
log1loglog)(1log),(
2
2
1
12
Θ=��
�
� +Θ=
��
�
�
�+−Θ=
��
�
�
�+−+Θ=
=
−
=
8
Falls r = 1 , so heißt A Tridiagonalmatrix.
Die Lösung eines Gleichungssystems Ax = b, in dem A eine Tridiagonalmatrix ist,
mittels der Gauß-Elimination bedeutet, dass in jedem Schritt nur ein Eliminationsfaktor
lk berechnet werden muss und nur ein Wert, für eine Zeile der Matrix, neu berechnet
werden muss. Der Eliminationsfaktor lk+1 kann aber erst dann berechnet werden,
nachdem die Zwischenmatrix A(k) bestimmt wurde, welche wiederum von lk abhängig
ist. Diese Datenabhängigkeit führt dazu, dass die Berechnung der oberen
Dreiecksmatrix sequentiell durchgeführt werden muss. Eine parallele Implementierung
der Gauß-Elimination ist bei Matrizen mit dieser Struktur also nicht sinnvoll.
Direkte Verfahren zur Lösung eines linearen Gleichungssystems mit einer solchen
Struktur sind z.B. rekursives Verdoppeln oder als eine Variante davon die zyklische
Reduktion. Beide Methoden sind für lineare Gleichungssysteme geeignet deren Matrix
A symmetrisch (d.h. A=AT) und positiv definit (d.h. xTAx > 0 für alle x∈Rn x≠0) ist,
außerdem darf der Vektor b nur von Null verschiedene Elemente enthalten.
Ausgangspunkt beider Verfahren ist das folgende Gleichungssystem:
Die Basis beider Verfahren ist es xi-1 und xi+1 aus der i-ten Gleichung, durch einsetzen
von geeigneten Vielfachen der Gleichungen i+1 und i-1 zu eliminieren. Dies führt zu
einen Gleichungssystem in dem die Nebendiagonalen von der Hauptdiagonalen
„weggeschoben“ werden. Nachdem diese Vorgehen zum ersten mal auf die
Tridiagonalmatrix angewandt wurde ergibt sich folgende Matrix:
nn
iiii
yb
ifürycba
ycb
=+==++
=+
+
x xa
1-n,2, x x x
xx
n 1-nn
1ii1-i
12111
�
��������
�
�
��������
�
�
=−
nn
n
ba
c
ba
cbcb
A
000
000
00
)1(
)1(2
3)1(
3
)1(2
)1(2
)1(1
)1(1
)1(
��
�
��������
�
�
��������
�
�
=
−
nn
n
ba
c
ba
cba
cb
A
0
0
1
33
222
11
��
� Abbildung 6 Tridiagonalmatrix
9
Wird dies nun entsprechend oft wiederholt, so entsteht eine Matrix auf der nur noch die
Elemente der Hauptdiagonalen ungleich Null sind und die Lösung abgelesen werden
kann. Zur Berechnung dieser Diagonalmatrix benötigt man bei einer n × n
Ausgangsmatrix �log n � Schritte. In jedem dieser Schritte k=1,..., �log n � werden die
Gleichungen i-2k-1, i, i+2k-1
betrachtet, um aus Gleichung i die Unbekannten xi+2^(k-1) und xi+2^(k-1) zu eliminieren.
Hierzu müssen jeweils die Werte ai(k) ,b i
(k),c i(k),y i
(k) für i = 1,...,n berechnet werden.
Da die Berechnung der αi(k), βi
(k) eine Division durch bi(k) erfordert sind beide
Verfahren, Rekursives Verdoppeln und zyklische Reduktion, nur dann Anwendbar,
wenn alle bi(k) ≠ 0. Für Insgesamt benötigt das Verfahren des rekursiven Verdoppelns
�log n� Schritte um die Lösung zu ermitteln, in jedem dieser Schritte werden O(n)
Werte berechnet, mit der daraus resultierenden Laufzeit O(n⋅log n) ist das Verfahren
eigentlich, gegenüber der Gauß-Elimination, im Nachteil. Diese benötigt für eine
Tridiagonalmatrix nur eine Laufzeit von O(n), da nur eine Eliminationsfaktor und ein
Koeffizient in jedem der n-1 Transformationsschritte berechnet werden muss. Der
Vorteil gegenüber der Gauß-Elimination liegt aber nun darin, dass die Berechnung der
Werte innerhalb eines Schrittes parallel ausgeführt werden kann, da die Berechnung von
ai(k) ,b i
(k),c i(k),y i
(k) nur von Werten aus Schritt k-1 abhängig ist . Rekursives Verdoppeln
)1(22
)1(22
)1(2
)1(2
)1(2
)1()1(2
)1(
)1(2
)1(22
)1(22
)1(2
11111
11
11111
−++
−++
−+
−+
−+
−−−
−
−−
−−−
−−−
−−
−−−−−
−−
−−−−−
=++
=++
=++
kii
kii
kii
ki
kii
kii
kii
ki
kii
kii
kii
ki
kkkkkk
kk
kkkkkk
yxcxbxa
yxcxbxa
yxcxbxa
2-n ,...,1ifür :
n, ,...,2ifür :
mit
, ,...,1ifür
, ,...,1ifür b
sonst, 0c und 2 ,...,1ifür
sonst, 0 und ,,...,12ifür
1-k)1(
2
)1()(
1-k)1(
2
)1()(
)1(2
)(1)-(ki
)1(2
)()(
)1(2
)(1)-(ki
)1(2
)()(
)(k)1(2
)()(
)(k)1(2
)()(
1
1
11
11
1
1
=−=
=−=
=⋅++⋅=
=⋅++⋅=
=−=⋅=
=+=⋅=
−+
−
−−
−
−+
−−
−+
−−
−+
−−
−
−
−−
−−
−
−
ki
kik
i
ki
kik
i
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
ki
k
k
kk
kk
k
k
b
c
b
a
nyyyy
naca
ncc
anaa
β
α
βα
βα
β
α
10
p 1,..., ifür
~~
~
~~
)(
)(
)(
)(
=
���
�
���
�
�
=====
⋅
⋅
⋅
⋅
⋅
qii
Qqii
Qqii
Qqii
Qqii
xx
yy
cc
bb
aa
und zyklische Reduktion unterscheiden darin, dass die zyklische Reduktion nicht alle
Zahlen ai(k) ,b i
(k),c i(k),y i
(k) in jedem Schritt berechnet, sondern die Zahl der zu
berechnenden Werte in jedem Schritt halbiert. Für k=1,...,�log n � werden die Zahlen
ai(k) ,b i
(k),c i(k), y i
(k) für i=2k,...,n mit Schrittweite 2k berechnet. Hierdurch wird aber
eine abschließende Substitutionsphase zur Berechnung der Lösung notwendig, siehe
Abb.7. In dieser werden für k= �log n �-1,...,0 die Elemente xi mit i=2k,...,n mit
Schrittweite 2k+1 durch:
berechnet.
Bei parallelen Realisierung der zyklischen Reduktion werden die Zeilen der n × n
Tridiagonalmatrix A blockweise auf p Prozessoren aufgeteilt. Zur Vereinfachung soll n
= p⋅q mit q∈ Ν und q = 2Q mit Q ∈ N sein. Jeder Prozessor Pi speichert nun einen
Zeilenblock der Größe q, der die Zeilen (i-1)q+1,...,i*q für 1 ≤ i ≤ p enthält, ab.
Der Algorithmus teilt sich in drei Phasen auf:
1. zyklische Reduktion: Prozessor Pi, mit 1 ≤ i ≤ p, berechnet die ersten Q =log q
Stufen der zyklischen Reduktion. Pi berechnet also in jeder Stufe k = 1, ...,Q
ai(k) ,b i
(k),c i(k),y i
(k) mit j=(i-1)*q + 2^k, ..., i*q mit Schrittweite 2k. Jeder Pi muss für
diese Berechnungen jeweils von Pi+1 (falls i<n) und Pi-1 (falls i>0) eine Zeile (i*q+1
bzw. (i-1)*q), die in der vorangegangenen Stufe berechnet wurde, empfangen.
2. rekursives Verdoppeln: Nach Abschluß der ersten Phase hat jeder Prozessor Pi die
Gleichung
mit p1,...,ifür ~~~~~~~11 ==++ +− iiiiiii yxcxbxa
)(2
)(2
)()(
ki
ik
iik
ik
ii
b
xcxayx
kk +− ⋅−⋅−= Gleichung 5
Abbildung 7 Schema der
zyklische Reduktion mit 4
Gleichungen nach Phase 1 liegt nur
x4 vor x2 wir im ersten Schritt der
Substitutionsphase ermittelt x1 und
x3 im zweiten
i=1
i=2
i=3
i=4
11
)()()()( ~,~,~
,~ ki
ki
ki
ki ycba
)'(
)'(~~
Ni
Ni
ib
yx =
)1()1()1()1( ~,~,~
,~ −−−− kj
kj
kj
kj ycba
vorliegen.Dieses neue nur noch p-dimensionale Gleichungssystem wird nun mittels
rekursiven Verdoppeln in �log p� Schritten gelöst. Alle Prozessoren haben dabei nur
noch eine Zeile zu bearbeiten. In jedem dieser k=1,..., �log p� Schritte berechnet
Pi die vier Koeffizienten , wofür er die Koeffizienten
aus dem Schritt k-1 für j∈{i-2k-1,i+2k-1} Pi-2^(k-1)
und Pi2^(k-1) benötigt. Abschließend errechnet Pi nach Schritt N=�log p�
3. Substitutionsphase: Nach dem jeder Prozessor Pi x(i⋅q) in Phase 2. berechnet hat,
werden in dieser Phase von jedem Prozessor Pi die noch fehlenden xj
mit j =(i-1)q+1,... ,(i⋅q)-1, entsprechend der Gleichung 5, berechnet. Diese Phase
erfordert Q Schritten, in jedem Schritt k=Q-1,...,0 berechnet die Prozessoren die
Elemente xj, j= 2k,...,n, und Schrittweite 2k+1 .Jeder einzelne Prozessor Pi berechnet
davon die xj für die gilt j div q +1 = 1.Hiefür benötigt er die Werte x(i-1)q, x(i+1)q von
sein Nachbarprozessoren Pi-1 und Pi+1 falls diese vorhanden sind. In Abb. 8 wird die
Vorgehensweise in den einzelnen Phasen noch einmal verdeutlicht.
Abbildung 8 Schematische Darstellung der zyklischen Reduktion mit n=8 p=2 q=4
Q=2, in der die Gleichungen durch Knoten symbolisiert werden.
i=1
i=2
i=3
i=4
i=5
i=6
i=7
i=8
12
In der ersten Phase mit Q= log q= log(n/p) Stufen beginnt jeder Prozessor mit q
Gleichungen und muss somit in der ersten Stufe maximal 4q/2 Werte berechnen. In
jeder weiteren Stufe halbiert sich die Anzahl an Gleichungen und somit an Werten die
jeder Prozessor berechnet. Für Stufe k=1,...,log q berechnet jeder Prozessor also 4q/2k
Zahlen, wozu 4 Operationen pro Zahl erforderlich sind. Daraus ergibt sich ein
Berechnungsaufwand von:
Darüber hinaus sendet und empfängt jeder Prozessor in jeder Stufe maximal 2
Nachrichten mit 4 Elementen. Der Kommunikationsaufwand hierfür beträgt:
In jeder Stufe der Phase 2, mit insgesamt �log p� Stufen, muss jeder Prozessor jeweils
pro Stufe eine Zeile mit 4 Werte berechnen (wiederum 4 Operationen pro Wert
notwendig). Dazu kommt noch eine Operation zur Berechnung von x(i⋅q) durch
Prozessor Pi. Der daraus resultierende Rechenaufwand beträgt:
Wie in Phase 1 empfängt jeder Prozessor 2 Nachrichten pro Stufe, mit jeweils 4 Werten.
Die dritte Phase benötigt weitere Q Stufen, in Stufe k=0,...,Q-1 berechnet jeder
Prozessor 2k Elemente des Lösungsvektors, wofür 5 Operationen notwendig sind.
Außerdem empfängt jeder Prozessor einen Wert von jedem seiner beiden
Nachbarprozessoren.
Insgesamt ergibt dies einen Berechnungsaufwand und Kommunikationsaufwand von:
)4(log2)4(2),( 221 ssss tpn
tQpnC ⋅⋅=⋅=
opopopop ttptptpnT +⋅��=+��⋅= log16log44),(2
)4(log2),( 22 sstppnC ⋅��=
���
�
�−⋅=−⋅=−⋅=⋅=
−
=
15)1(5)12(525),(1
03 p
ntqtttpnT opop
Qop
Q
k
kop
)1(2),( 23 sstpnC ⋅=
op
Q
kkop t
pnq
tpnT ⋅≤⋅= =
1624
4),(1
1
13
Damit verteilt sich der Berechnungsaufwand fast (bis auf einen kleinen Zusatz von
16⋅log p) optimal auf die einzelnen Prozessoren und für große n wird der linear
ansteigende Berechnungsaufwand den nur logarithmisch ansteigenden
Kommunikationsaufwand übersteigen.
3 Iterative Verfahren
3.1 Klassische iterative Verfahren
Für dünnbesetzte Matrizen ist das Gauß-Elimination nur bedingt geeignet, da die
Umformung der Matrix A zur Auffüllung mit Nichtnullelementen führen kann .
Entweder muss ein Algorithmus verwendet werden, der an die spezielle Matrixstruktur
angepaßt ist, z.B. zyklische Reduktion oder man verwendet ein iteratives Verfahren.
Diese liefern, im Gegensatz zu den direkten Verfahren keine exakte Lösung, sondern
berechnen vielmehr eine Näherung an die gesuchte Lösung. Für das Gleichungssystem
Ax =b mit Matrix A∈ Rn×n und Vektor b∈ Rn wird eine Folge von Approximations-
vektoren {x(k)}k=1,2,..., berechnet die gegen die gesuchte Lösung x*∈Rn konvertieren. Ein
wichtiges Kriterium zur Bewertung von Iterationsverfahren ist daher die Konvergenz-
geschwindigkeit. Basis klassischer Iterationsverfahren ist die Zerlegung der Matrix A in
A=M-N mit M,N∈ Rn×n . M ist dabei eine nichtsinguläre Matrix und M-1 sollte
möglichst leicht zu berechnen sein ,z.B. eine Diagonalmatrix. Die gesuchte Lösung x*
des Gleichungssystems erfüllt die Gleichung Mx* = Nx* +b. Daraus ergibt sich auch die
Iterationsvorschrift Mx(k+1) = Nx(k) + b. Üblicherweise wird diese aber in folgender
Schreibweise angegeben:
mit C:=M-1 N und d:=M-1 b. Das Iterationsverfahren heißt konvergent gegen x* ,wenn
ein x* existiert, so daß {x(k)}k=1,2,..., unabhängig von der Wahl des Startvektors x(0) ∈ Rn
gegen x* konvergiert.
dCxx kk +=+ )()1( Gleichung 6
� �
� � ).1(22)4(2log2)1(22)4(2log2log2),(
log162145log1616),(
stsstsnstsstsppn
pnC
tppn
tpn
ppn
pnT opop
⋅+⋅⋅≅⋅+⋅���
�
�+⋅=
⋅���
�
�⋅+≅⋅��
�
�
�−+⋅+=
14
3.1.1 Jacobi-Verfahren
Beim Jacobi-Verfahren erfolgt die Zerlegung der Matrix A in A=D-L-R (D,L,R∈ Rn×n),
hierbei entspricht D der Diagonalmatrix, L der unteren Dreiecksmatrix und R der
oberen Dreiecksmatrix, jeweils ohne Diagonale. Die Zerlegung wird in Form der
Iterationsvorschrift Dx(k+1) =(L+R)x(k) + b genutzt. Mit CG:=D-1(L+R) bzw.
ergibt sich der Iterationsschritt in Komponentenschreibweise:
Um zu überprüfen ob das Jacobi-Verfahren konvergiert wird häufig als Kriterium das
einfach zu überprüfende starke Zeilensummenkriterium:
verwendet. Erfüllt Matrix A dieses Kriterium, dann konvertiert das Jacobi-Verfahren.
Das Verfahren bricht ab wenn die Abweichung von x(k) zur Gesuchten Lösung x(*)
kleiner ist, wie die festgelegte Toleranzgrenze ε. Da x(*) natürlich nicht bekannt ist wird
der relative Fehler als Abbruchkriterium genutzt.
wobei ||.|| eine frei zu wählende Vektornorm bezeichnet z.B. ||x||∞ = maxi=1,...,n|xi| oder
||x||2=( n i=1|x|2)½ .
Eine Möglichkeit der parallelen Implementierung des Jacobi-Verfahrens ist die
Verwendung von Parallelisierungen der Matrix-Operationen, die für die Berechnung
von x(k+1) benötigt werden. Dies sind eine Matrix-Vektor-Multiplikation von (L+R) mit
x(k) , eine Vektor-Vektor-Multiplikation des Ergebnisses aus 1. mit Vektor b sowie, eine
Matrix-Vektor-Multiplikation mit D-1.
Die andere Möglichkeit der parallelen Implementierung besteht darin, dass die
Berechnung der einzelnen Elemente des Vektors xi(k+1),für einen Iterationsschritt
unabhängig von einander ist. Dies ermöglicht der Verwendung von maximal p=n
Prozessoren. Bei Verwendung von p Prozessoren speichert jeder Prozessor Pi mit
i=1,...,p die Werte xi(n/p),..., x(i+1)(n/p)-1 inklusive der dazugehörigen Zeilen der Matrix A
��� ≠−
== = sonst. 0, ijfür /
cmit )( ij,...,1,iiij
njiijGaa
cC
n, 1,...i , 1
,1
)()1( =��
�
�
�−=
≠=
+n
ijj
kjiji
ii
ki xab
ax
n , 1,...i, ,1
=> ≠=
n
ijjii aija
)1()()1( ++ ≤− kkk xxx ε
Gleichung 7
15
und des Vektor b. Jeder Prozessor Pi führen nun alternierend folgende Aktionen aus bis
das Abbruchkriterium erfüllt ist:
1. Da x(k) zur Berechnung der neuen Approximation x(k+1) notwendig ist, sendet
Prozessor Pi seine lokal gehaltenen Elemente des Vektors x(k) an alle anderen
Prozessoren, dies erfordert einen All-to-All Broadcast.
2. Prozessor Pi berechnet gemäß Formel 7 die neue Approximation x(k+1) für seine
Elemente xj(k+1) mit j= i(n/p),...,(i+1)(n/p)-1
3. Prozessor Pi berechnet lokal, für das Abbruchkriterium, die Differenzen xj(k+1) - xj
(k)
für alle j= i(n/p),...,(i+1)(n/p)-1. Diese lokalen Ergebnisse müssen dann zusammen
geführt werden, um den Test entsprechend der gewählten Vektornorm
durchzuführen.
In Schritt 1. sendet Pi n/p Werte an p-1 Prozessoren dies ergibt einen
Kommunikationsaufwand von �((p-1) ⋅ (n/p)). In Schritt 2. muss jeder Prozessor n/p
Elemente des Vektors x(k+1) berechnen, zur Berechnung eines Elementes muss gemäß
Formel 7 eine Summe von j=1,...,n gebildet werden, damit ergibt sich ein
Rechenaufwand von �(n⋅ (n/p)) für Schritt 2. Der Test auf Abbruch in 3. erfordert einen
Aufwand von � (log p). Nachteil des Verfahrens ist , dass wenn es konvergiert es eine
niedrige Konvergenzrate hat. Wieviel Iterationen genau benötigt werden bis das
Verfahren konvergiert ist für den Allgemeinen Fall schwierig vorherzusagen, es wurde
aber gezeigt, dass e(n2/2) Iterationen notwendig sind damit der Fehler um 10-e sinkt.
Für n= 64 wären ca. 3(642/2)= 6144 Iterationen notwendig um den Fehler um 10-3 zu
reduzieren. Die Konvergenzgeschwindigkeit läßt sich aber durch einführen eines
sogenannten Relaxtionsparameters ω ∈ R verbessern. Der Iterationsschritt in
Komponentenschreibweise lautet dann wie folgt
Diese Variante des Jacobi-Verfahrens wird auch als JOR–Verfahren bezeichnet. Die
optimale Wahl des Parameters ist aber von den numerischen Eigenschaften des zu
lösenden Gleichungssystems abhängig und die Bestimmung daher nicht immer einfach.
n, 1,...i , )-(1 )(
,1
)()1( =+��
�
�
�−=
≠=
+ ki
n
ijj
kjiji
ii
ki xxab
ax ωω
16
3.1.2 Gauß-Seidel-Verfahren
Das Gauß-Seidel-Verfahren (GS-Verfahren) nutzt die gleiche Zerlegung der Matrix A,
wie das Jacobi-Verfahren A=D-L-R (D,L,R∈ Rn×n) auch Konvergenzkriterium und
Abbruchkriterium sind die Gleichen des Jacobi-Verfahrens. Die Zerlegung wird aber in
anderer Form für den Iterationsschritt genutzt (D-L)x(k+1)=Rx(k)+b. Iterationsmatrix CE
ist somit definiert als CE:=(D-L)-1R und der Iterationsschritt lautet in
Komponentenschreibweise
Durch Einbeziehung der neu Berechneten Approximation xj(k+1) , j=1,...,i-1 zur
Berechnung von xi(k+1) i=1,...,n wird eine, gegenüber dem Jacobi-Verfahren, deutlich
verbesserte Konvergenzrate erzielt. Diese läßt sich, ebenso wie beim Jacobi-Verfahren
noch einmal verbessern durch die einführen eines Relaxtionsparameters ω ∈ R. Der
Iterationsschritt in Komponentenschreibweise ergibt sich wie folgt:
In dieser auch als SOR-Verfahren bezeichneten Variante wird der Relaxtionsparameter
zur Kombination von altem und neuem Iterationsvektor genutzt. Die Konvergenzrate
des GS-Verfahren liegt bei n⋅(e/3) Iterationen um den Fehler um den Faktor 10-e zu
senken, für n=64 wären ca. (64)3/3=64 Iterationen notwendig um den Fehler um 10-3 zu
reduzieren, gegenüber 6144 Iterationen im Jacobi-Verfahren. Allerdings führt diese
Einbeziehung der neu berechneten Approximation xj(k+1) j=1,...,i-1 zur Berechnung von
xi(k+1) i=1,...,n zu Datenabhängigkeiten die eine Parallelisierung des Verfahrens
erschweren. Die Berechnung von der xi(k+1) i=1,...,n erfordert zuvor die Bestimmung
aller xj(k+1) für j=1,...,i-1, dadurch kann die Berechnung der xi
(k+1) i=1,...,n nicht parallel
erfolgen. Parallelität ist nur innerhalb der Berechnung eines xi(k+1) i=1,...,n, durch
Berechnung von Teilsummen auf verscheiden Prozessoren möglich. Für dichtbesetzte
Matrizen ist das Verfahren daher ungeeignet. In Dünnbesetzte Matrizen A=(aij)ij=1,...,n ,
mit vielen aij=0 ergeben sich weniger Datenabhänigkeiten (ist aij=0 dann kann xi(k+1)
ohne Kenntnis von xj(k+1) berechnet werden, auch wenn j<i). Die Datenabhängigkeiten
lassen sich durch in einem Präzendenzgraphen darstellen, gerichtete Kante zwischen
Knoten j und Knoten i, falls j zur Berechnung von i benötigt wird(also aij ≠0). Abb.9
n,1,i, 1 1
1
)(
1
)1()1(�=
��
�
�
�−−=
−
= +=
++i
j
kj
n
ijij
kjiji
ii
ki xaxab
ax Gleichung 8
n,1,i, )x-(1 (k)i
1
1
)(
1
)1()1(�=+
��
�
�
�−−=
−
= +=
++ ωω i
j
kj
n
ijij
kjiji
ii
ki xaxab
ax
17
links zeigt einen resultierenden Präzendenzgraphen der ersten zwei Schritte des GS-
Verfahrens. Hieraus ist zu erkennen, das nur x3 und x4 parallel berechnet werden
können.
Die Parallelität läßt sich aber durch Umordnung der Gleichungssystems erhöhen. Hierzu
wird im Beispiel aus Abb.9 rechts die Numerierung der Variablen geändert, 2 → 4, 4 →
3, und 3 → 4. Aus dem daraus resultierenden Präzendenzgraph (Abb. 9 rechts) ist zu
erkennen, dass nun eine parallel Berechnung von x1, x2 und x3 möglich ist. Damit stellt
sich für das GS-Verfahren die Frage, nach einer guten Anordnung bzgl. größt möglicher
Parallelisierung. Eine Möglichkeit ist die sogenannte Rot-Schwarz-Anordnung, die sich
für Gleichungssysteme eignet die von einer Gitterstruktur herrühren. Solche
Gitterstrukturen treten z.B. bei der numerischen Lösung partieller Differential-
gleichungen auf oder wie im Beispiele aus Abbildung 10 bei der Bestimmung des
Temperaturverlaufs in einem Wasserbad.
Die Temperatur in Punkt xij wird dabei mit Hilfe der vier Nachbartemperaturen
bestimmt xi,j= (xi-1,j + xi+1,j + xi,j-i + x i,j+1) / 4. Für n = 16 Punkte ergibt sich folgendes
Gitter wie in Abbildung 11 zu sehen, die Struktur der Matrix A, des daraus
resultierenden Gleichungssystems ist ebenfalls in Abb. 11 dargestellt.
x1 x2 x3 x4 x1 x2 x3 x4Abbildung 9Beispiel für einen
resultierenden
Präzendenzgraphen
des GS-Verfahrens
Abbildung 10 Vorgabe für
die Bestimmung des
Temperaturverlaufs in einem
Wasserbad
18
��
�
�
�=��
�
�
����
�
�=⋅
���
�
�=��
�
�
�=��
�
�
�=
2
1
S
R
S
R
S
R
ˆˆ
x
ˆ.
DEFD
ˆA
00F-0
R 0E-00
L ,D00D
D
b
bxx
Diese Gleichungssystem Ax=b mit A 16×16 Matrix ist für die parallele
Implementierung des GS-Verfahren nicht brauchbar, z.B. neu Approximation von x6
abhängig von Berechnung der neuen Approximation von x2 und x5. Um eine
geeignetere Anordnung zu erhalten werden die Punkte in Rote und Schwarze Punkte
aufgeteilt. Die Roten werden von 1,...,nS, die Schwarzen von nS+1,... ,nR+nS
durchnumeriert (nR+nS=n). Durch umsetzen der Forderung, dass rot Punkt nur
schwarze
Nachbarn haben dürfen und umgekehrt ergibt sich ein Gitter wie in Abb. 12. Die
Struktur der Matrix des aus diesem Gitter abgeleiteten Gleichungssystems ist ebenfalls
in Abb. 12 zusehen. Diese neue Matrix Â, die eine Permutation der Ausgangsmatrix ist,
besteht aus den zwei Diagonalmatrizen DR, DS und den dünnbesetzten Bandmatrizen E,
F. Die Zerlegung der Matrix  erfolgt gemäß GS-Verfahren in  =D-L-R (D,L,R∈
Rn×n) mit
�����������������������
�
�
�����������������������
�
�
xxxx
xxxxxxx
xxxxxxxxx
xxxxxxx
xxxxxxx
xxxxxxxxx
xxxxxxxx
xxxxxxxxx
xxx
Abbildung 11 Beispiel
eines 4*4 Gitters und der
schematischen Darstellung
der dazugehörigen Matrix
Abbildung 12 Gitter
in Rot – Schwarz -
Anordnung und dazu
gehörige Matrix A
�����������������������
�
�
�����������������������
�
�
xxxxxxx
xxxxxx
xxxxxxxxx
xxxxxxxxx
xxxxxxxxx
xxxxxxxxx
xxxxxxx
xxxxxxx
19
Dies führt zu dem Iterationsschritt
oder in Komponentenschreibweise, wobei N(i) die Menge der vier Nachbargitterpunkte
von Gitterpunkt i bezeichnet:
Jetzt wird deutlich ,dass xR(k+1) nur noch von xS
(k) abhängt und xS(k+1) von xR
(k+1).
Abwechselnd können nun die roten Iterationsvektoren xR(k+1) und die schwarzen
Iterationsvektoren xS(k+1), jeweils parallel, berechnet werden. Damit wird für die roten
Iterationsvektoren, bzw. die schwarzen Iterationsvektoren, jeweils der gleiche
Parallelitätsgrad, wie beim Jacobi-Verfahren erreicht und dies mit der größeren
Konvergenzgeschwindigkeit des GS-Verfahrens. Für die Parallele Implementierung
können maximal p=nr (bzw. p=ns) Prozessoren verwendet werden. Die Aufteilung der
zur Berechnung der Approximationen xi wird anhand des Gitters vorgenommen. Wird
Prozessor pj Gitterpunkt i zugeordnet so ist er für die Berechnung der Approximationen
von xi zuständig. Die Aufteilung der Gitterpunkte kann z.B. Zeilen blockweise erfolgen,
bei p Prozessoren erhält jeder Prozessor √n / p Gitterzeilen und damit ½(n / p) rote und
schwarze Gitterpunkte. Ausgehen von einem Startvektor berechnen die Prozessoren
folgende Schritte, die solange iteriert werden bis der Test auf Abbruch erfolgreich ist..
1. Da Vektor xS(k) zur Berechnung der neuen Approximation von xR
(k+1) notwendig ist,
sendet Prozessor Pi seine lokal gehaltenen Elemente des Vektors xS(k) an alle
anderen Prozessoren, dies erfordert einen All-to-All Broadcast.
2. Prozessor Pi berechnet gemäß Formel 9 die neue Approximation xR(k+1) für seine
Elemente von xR(k)
3. Prozessor Pi sendet seine lokal gehaltenen Elemente des Vektors xR(k+1) an alle
anderen Prozessoren (All-to-All Broadcast).
( ) ( )
( ) ( ) R)(
)1(,
,
)1(
S)(
)()1(
n, 1,...i , ˆˆ
1
,n, 1,...i , ˆˆ1
=����
�
�
�
−=
=��
�
�
�−=
∉
+++
++
+
∉
+
n
iNjj
kRjnini
ninii
kS
iNjj
kSiji
iii
kR
xaba
x
xaba
x
ss
ss
Gleichung 9
1)(kR2
1)(kSS
(k)S1
1)(kRR
ˆD
FˆD++
+
⋅−=⋅
⋅−=⋅
xEbx
xbx
20
4. Prozessor Pi berechnet gemäß Formel 9 die neue Approximation xS(k+1) für seine
Elemente von xS(k)
5. Prozessor Pi berechnet lokal, für das Abbruchkriterium, die Differenz aus der alten
und neuen Approximation. Diese lokalen Ergebnisse müssen dann zusammen
geführt werden um den Test durchzuführen.
Der Rechenaufwand ist damit fast identisch zu dem des Jacobi-Verfahrens . In Schritt 2
und 4 muss jeder Prozessor jeweils n/2p Elemente der neuen Approximation berechnen,
was insgesamt die gleiche Anzahl, wie beim Jacobi-Verfahrens ergibt. Ein Unterschied
zum Jacobi-Verfahren besteht nur darin, dass eine All-to-All Broadcastoperation mehr,
Schritt 1 und 3, erforderlich ist . Dieser zusätzliche Aufwand wird aber durch dass
verbesserte Konvergenzverhalten mehr als kompensiert.
3.2 Methode der konjugierten Gradienten
Das Verfahren der konjugierten Gradienten (CG-Verfahren) ist eine iteratives Verfahren
zur Lösung von Gleichungssystemen der Form Ax =b mit Matrix A∈Rn×n positiv definit
(d.h. xTAx > 0 für alle x∈Rn x≠0) und symmetrisch(d.h. A=AT). Im Gegensatz zu
klassischen Iterationsverfahren, bei denen die Anzahl der Iterationsschritte vorher nicht
bekannt ist, liegt beim CG-Verfahren, bei exakter Rechnung, die Lösung nach n
Schritten vor, Rundungsfehler führen allerdings dazu, dass weitere Iterationen zur
Bestimmung einer hinreichend genauen Lösung benötigt werden. Für dünnbesetzte
Matrizen gilt dies nicht hier wurde gezeigt, dass eine hinreichend genaue Lösung in
wenigen als n Schritten berechnet werden kann. Die Lösung x* des oben beschriebenen
Gleichungssystems entspricht dem Minimum der Funktion
f(x) ist konvex und besitzt eindeutiges Minimum x*. Zur Bestimmung der Lösung x*
werden ausgehend von einem frei zu wählenden Startvektor x(0) weitere
Approximationen x(k) des Minimums ermittelt bis x(k) nahe genug an der gesuchten
Lösung x*. Im Iterationsschritt wird x(k+1) wie folgt bestimmt :
)()()1( )( kkk pkxx α+=+
xbAxxxf TT −=21
)(
21
Parameter α(k) entspricht dabei der Schrittweite und p(k) der Richtungsänderung. Die
Bestimmung von α(k) bei bekanntem p(k) erfolgt so, dass
Die Notwendige Bedingung für das Minimum lautet (im folgenden ohne Index k):
mit Residualvektor r =Ax-b. Daraus Folgt vermeide die Wahl von p, so dass
pT(Ax-b)=0. Die Hinreichende Bedingung für das Minimum
ist immer erfüllt da A, wie oben gefordert positiv definit ist. Die Unterschiedliche Wahl
von p(k) führt zu verschiedenen Iterationsverfahren, einfachste Möglichkeit ist Methode
des steilsten Anstieges mit p(k)= -r(k) (r =Ax(k)-b), diese hat aber eine sehr schlechte
Konvergenz. Besser geeignet ist das Verfahren der konjugierten Gradienten die
Richtungsvektoren p(k) sind dabei so zu wählen, dass sie entsprechend dem gestellten
Problem orthogonal zu einander sind, derartige Vektoren werden als A - orthogonal
oder konjugiert bezeichnet. Zwei Vektoren p,q ∈Rn sind bzgl. einer symmetrischen,
nicht singulären Matrix A konjugiert falls gilt qTAp=0. Ist A positiv definit so bilden die
linear unabhängigen Vektoren p1,...,pn, mit pi ≠0, i=1,...,n und (pi)TApj=0 für alle i,j eine
konjugierte Basis des Rn bzgl. A. Falls die Richtungsvektoren zur Minimierung eine
konjugierte Basis des Rn bzgl. A, so führt das oben beschriebene Iterationsverfahren zur
Minimierung bei exakter Rechnung nach maximal n Schritten zur exakten Lösung x*.
Der Algorithmus CG-Verfahrens sieht nun im einzelnen wie folgt aus:
Initialisierung: wähle beliebige Belegung von Startvektor x(0) und setzte
p(0) =-r(0) =b-Ax(0) sowie k=0
Iteration: solange ||r(k)||>ε berechne
1. Matrix - Vektor- Produkt: q(k) =Ap(k)
2. Schrittweite bestimmen: αk = [(r(k))T r(k)] / [ (p(k))T q(k) ]
3. Neue Näherung an Lösung ermitteln: x(k+1) =x(k) + αkp(k)
4. Update des Residualvektor: r(k+1) =r(k) + αkq(k)
))((min))(( )()( += + ∈ kk
Rk pxfpkxf αα α
App
rp
App
bAxpbAxpApp
dapxdf
T
T
T
TTT −= −−=⇔=−+=+ )(
0)()(
minααα
>=+ 0
)(2
2
Appda
pxdf Tα
22
5. Korrekturfaktor bestimmen βk= [(r(k+1))T r(k+1)]/ [(r(k))T r(k)]
6. Neue Suchrichtung berechnen: p(k+1) =-r(k+1)+ βkp(k)
7. Index erhöhen: k++
Die Parallele Implementierung des CG-Verfahren ergibt sich aus der Parallelisierung
der im Algorithmus verwendeten Basisoperationen. Diese sind im einzelnen eine
Matrix-Vektor-Multiplikation in 1., zwei innere Produkte in Schritt 2. (eins [(r(k))T r(k)]
wurde aber schon in der vorherigen Iteration in Schritt 5. berechnet (� γk=(r(k))T r(k)
zwischenspeichern), zwei Vektor-Additionen in Schritt 3. und 4., zwei innere Produkte
in Schritt 5. eins davon aus vorheriger Iteration bekannt und eine Vektor-Additionen in
Schritt 6. Insgesamt werden eine Matrix-Vektor-Multiplikation, zwei innere Produkte
und drei Vektor-Additionen in jedem Iterationsschritt berechnet. Bei der
Parallelisierung dieser Basisoperationen ist es wichtig, dass die Datenverteilung der
Matrix A und der Vektoren „zusammenpassen“, da sonst unnötige Kommunikation
erforderlich wird. Eine Möglichkeit der Realisierung ist die zeilenblockweise
Verteilung der Matrix A und blockweise die Verteilung der Vektorelemente der
Vektoren x(k) , p(k) , r(k) , q(k) .Der Algorithmus für jeden Iterationsschritt k laut wie
folgt:
1. All-to-All Broadcast zur Verteilung p(k)
2. parallele Matrix – Vektor – Multiplikation zur Berechnung von q(k)
3. parallele Berechnung des inneren Produktes (p(k))T q(k) und von αk
4. lokale Berechnung von x(k+1) =x(k)+ αkp(k) und r(k+1)=r(k)+ αkq(k)
5. parallele Berechnung des inneren Produktes (r(k+1))T r(k+1) und βk
6. lokale Berechnung von p(k+1)=-r(k-+)+ βkp(k)
Die Laufzeit eines Iterationsschrittes des Algorithmus ,bei p Prozessoren, setzt sich
zusammen aus der Laufzeit zur Berechnung von drei Additionen oder auch axpy-
Operation genannt (a x plus y ). Die Parallele Implementierung einer axpy-Operation
erfordert einen Aufwand von
���
�
�Θ=⋅⋅=
pn
tpn
T opsaxpy 2
23
sowie dem Berechnungsaufwand für zwei innere Produkte, die Prozessoren bilden
jeweils die lokalen Summen, die dann über eine Reduktion zusammen geführt werden
mit einem Aufwand der sich zu:
Abschätzen läßt. Und der Laufzeit einer Matrix-Vektor-Muliplikation dafür muß
zunächst p(k) an alle Prozessoren verteilt werden, dann Berechnung der Skalarprodukte
ajT p(k) lokal von jedem Prozessor für alle lokalen Zeilen aj die er besitzt. Der Aufwand
beträgt hierfür:
Insgesamt ergibt dies einen Aufwand für das CG – Verfahren in jedem Iterationsschritt
von:
Da bei rundungsfehlerfreier Rechnung nach höchstens n Schritten mit der exakten
Lösung vorliegt und für dünnbesetzte Matrizen schon in weniger ,als n Schritten eine
hinreichend genaue Lösung ermittelt werden kann ist das Verfahren anderen
vorzuziehen wenn die Voraussetzungen, Matrix A symmetrisch und positiv definit,
erfüllt sind.
4 Zusammenfassung
Das Gauß-Eliminations-Verfahren eignet sich nur für fast voll besetze Matrizen. Seine
Vorteile liegen in der Berechnung einer exakte Lösung in Vorhersagbarer Rechenzeit
außerdem ist das Verfahren auf fast alle linearen Gleichungssystem Anwendbar. Für
dünnbesetzte Matrizen ist das Gauß-Eliminations-Verfahren, aufgrund des entstehenden
fill-in, hingegen wenig geeignet. In diesem fall sind andere Verfahren vorzuziehen
wenn deren Voraussetzungen erfüllt sind. Linearen Gleichungssystem mit
Tridiagonalmatrizen können durch zyklische Reduktion gelöst werden, als direktes
���
�
�+Θ= p
pn
T ip log
)()(2
pn
nnpn
ppn
T MVM +Θ=+Θ=
( )
��
�
�
�+++Θ=
��
�
�
����
�
�+��
�
�
�+++Θ=⋅+⋅+Θ=
pnpn
pn
pn
ppn
pn
nTTTT saxpyipMVMCG
log5
3log232
2
2
24
Verfahren liefer es eine exakte Lösung. Ansonsten sind Iterative Verfahren
anzuwenden, wobei das GS-Verfahren dem Jacobi-Verfahren vorzuziehen ist, da das
schneller konvergiert, insbesondere in der Variante SOR. Für Gleichungssysteme mit
positiv definiter Koeffizientenmatrix ist die Methode der konjugierten Gradienten
anderen Verfahren vorzuziehen, im Gegensatz zu den klassischen Iterationsverfahren,
die zur exakten Berechnung des Lösungsvektors eine nicht vorhersagbare Anzahl von
Iterationen berechnen müssen, stoppt das CG-Verfahren bei rundungsfehlerfreier
Rechnung nach höchstens n Schritten mit der exakten Lösung.
25
A Anhang Implementierung der Gauß-Elimination
Implementierung der zeilenorientieren Gauß-Elimination in Pseudocode für den
Hyperwürfel:
local n {Dimension des Gleichungssystems}a[n/p][n] {Koeffizienten der Matrix}b[n/p] {rechte Seite des Gleichungssystems}marked[n/p] {markieren der Pivotzeile}pivot[n/p] {Iteration in der Zeile Pivotzeile war}pickedmagintude {Pivot value}winner {Id des Prozessor der Pivotzeile besitzt}k,i,r
for all Pid, where 1 ≤ id ≤ p {
// marked initialisierenfor (k = 0; k< n/p; k++) {marked[i]=0;}
for (k = 0; k< n-1; k++) {
//lokales Maximum für Pivotzeile ermittelnmagnitude = 0;for(i=0; i< n/p;i++) {
if(marked[i]=0 && abs(a[i][k])> magnitude){picked=i; magnitude= abs(a[i][k]);
}}//globales Maximum für Pivotzeile ermittelnwinner = MAX_TOURNAMENT(id,magnitude,id);
// Gewinner markiert Pivotzeile ...if (winner == id){
marked[picked]= 1; pivot[picked]=k;for(i=k+1;i<n;i++){
tmp_vektor [i]=a[picked][i];}tmp_vektor[n]=b[picked]
}// ... und verteilt sie an die AnderenHYPERCUBE.BROADCAST(id,winner,
tmp_vektor[k+1...n])
// Neu Berrechnung der Matrix-Elementefor(i = 0; i < n/p; i++){
if(marked[i]=0){tmp=a[i][k]/tmp_vektor[k];for (r = k+1;r<n;r++){
a[i][r]= a[i][r]-tmp_vektor[r]*tmp}
}}
26
// Für die Zeile die nie Pivotzeile war ...for(i=0; i< n/p;i++){
// ...pivot setztenif(marked[i]=0){pivot[i]=n-1;}
}}
27
Literaturverzeichnis
[Qu94] Quinn, M.J.: Parallel Computing Theory and Practise, McGraw-Hill 1994.
[BU01] Buchholz , P.: Skript zur Vorlesung Verteilte numerische Algorithmen im
SS01, Abrufbar unter : http://www.iai.inf.tu-dresden.de/ms/
[KU00] Kuchen H. Skript zur Vorlesung Parallele Algorithmen im WS 00
[Ra98] Rauber, T. und G. Rünger: Parallele und verteilte Programierung,
Springer 1998.
top related