Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer...

31
Westfälische Wilhelms-Universität Münster Ausarbeitung Lösung linearer Gleichungssysteme (im Rahmen des Seminar „Parallele Programmierung“) Nico Ziborius Themensteller: Prof. Dr. Herbert Kuchen Betreuer: (Prof. Dr. Herbert Kuchen) Institut für Wirtschaftsinformatik Praktische Informatik in der Wirtschaft

Transcript of Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer...

Page 1: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 2: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

Westfälische Wilhelms-Universität Münster

Page 3: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 4: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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.

Page 5: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 6: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 7: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 8: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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.

Page 9: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 10: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Θ=��

� +Θ=

��

�+−Θ=

��

�+−+Θ=

=

=

Page 11: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 12: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

β

α

βα

βα

β

α

Page 13: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 14: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 15: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 16: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

⋅+⋅⋅≅⋅+⋅���

�+⋅=

⋅���

�⋅+≅⋅��

�−+⋅+=

Page 17: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 18: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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 ωω

Page 19: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 20: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 21: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 22: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 23: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

)(

Page 24: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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α

Page 25: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 26: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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

Page 27: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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.

Page 28: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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}

}}

Page 29: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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;}

}}

Page 30: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene

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.

Page 31: Lösung linearer Gleichungssysteme · PDF file1 1 Einleitung Die Lösung linearer Gleichungssysteme bildet den Kern vieler numerischer Anwendungen. Viele naturwissenschaftliche Phänomene