1 / 53

Lösung linearer Gleichungssysteme

Lösung linearer Gleichungssysteme. Seminar “Parallele Programmierung“ Nico Ziborius. Gliederung. 1. Einleitung. 2. Direkte Verfahren. 2.1 Gauß-Elimination. 2.2 zyklische Reduktion. 3. Iterative Verfahren. 3.1 klassische iterative Verfahren. 3.1.1 Jacobi-Verfahren.

howie
Download Presentation

Lösung linearer Gleichungssysteme

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Lösung linearer Gleichungssysteme Seminar “Parallele Programmierung“ Nico Ziborius

  2. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

  3. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

  4. 1. Einleitung • linearer Gleichungssysteme Kern vieler numerischer Anwendungen (z.B numerische Lösung partieller Differentialgleichungen) • Gesucht Lösung für Ax =b • nichtsinguläre Matrix ARnn • Vektor bRn und Lösungsvektor xRn • Direkte-Verfahren • exakte Lösung • z.B. Gauß-Elimination, zyklische Reduktion • Iterative-Verfahren • Näherung an die Lösung • z.B. Jacobi-Verfahren, Gauß-Seidel-Verfahren, Methode der konjugierten Gradienten

  5. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

  6. 2.1 Gauß-Elimination Idee: • Transformiere Ax=b in äquivalentes Gleichungssystem A‘x = b‘, so dass A‘ eine obere Dreiecksmatrix bildet.

  7. 2.1 Gauß-Elimination • Transformation benötigt n-1 Schritte • Schritt k , 1 k n-1 • A(k+1), b(k+1) aus A(k), b (k) berechnen: • Eliminationsfaktoren berechnen: • Elemente von A und b für k< j  n und k <i  n neu berechnen:

  8. 2.1 Gauß-Elimination • Lösung durch Rückwärtseinsetzen, für k = n,n-1,...,1: • sequentielle Implementierung Θ(n3) • (drei ineinander verschachtelte Schleifen)

  9. 2.1 Gauß-Elimination • Elemente auf der Diagonalen gleich Null • Fehler, da Division durch Null • sehr kleine Elemente auf der Diagonalen • Rundungsfehler • daher Elemente auf der Diagonalen durch ein andere, sogenannte Pivotelemente, ersetzten: • Spaltenpivotsuche: • betragsgrößte ark(k) aus akk(k),..,ank(k) suchen • vertauschen der Zeilen k und r , falls r  k.

  10. 2.1 Gauß-Elimination parallelen Implementierung: (Hypercube) • zeilenblockweise Aufteilung der Matrix A und des Vektor b • mit p = Anzahl Prozessoren, wird Pi für 1  i  p, ein n/p großer Zeilenblock zugewiesen • Jeder Prozessor Pi hat • ein 2-dimensionales Array a[n/p][n] (Koeffizienten Matrix A) • Array b[n/p] (Werten des Vektors b) • Array marked[n/p] (markieren der Pivotzeile) • Array pivot[n/p], (wann war Zeile Pivotzeile)

  11. marked pivot a b b 1 1 2 2 3 3 -4 -4 -6 -6 1 1 1 1 P1 P1 0 0 4 4 4 4 -1 -1 14 14 0 0 - - 0 0 -8 -8 -2 -2 10 10 3 3 0 0 - - P2 P2 0 0 -1 -1 -5 -5 11 11 7 7 0 0 - - 2.1 Gauß-Elimination • zunächst Array marked mit 0 initialisieren • dann in n-1 Schritten die obere Dreiecksmatrix berechnen: 1. lokales Pivotelement bestimmen: Pi ermittelt unter den Elementen a[r][k], mit r= 1...n/p und marked[r]=0 das lokale Maximum

  12. 2.1 Gauß-Elimination 2. globales Pivotelement bestimmen:Aus den lokalen Maxima wird das globale Maximum, durch Auruf der Funktion MAX_TOURNAMENT(int i, int value, int winner) ermittelt • 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;}

  13. 1 2 3 -4 -6 1 1 P1 0 4 4 -1 14 0 - 0 -8 -2 10 3 1 2 1 1 1 2 2 2 3 3 3 -4 -4 -4 -6 -6 -6 P2 P1 P1 P1 0 -1 -5 11 7 0 - 1 2 3 -4 -6 1 1 0 0 0 4 4 4 4 4 4 -1 -1 -1 14 14 14 P1 0 4 4 -1 14 0 - 0 0 0 -8 -8 -8 -2 -2 -2 10 10 10 3 3 3 P2 P2 P2 0 -8 -2 10 3 1 2 0 0 0 -1 -1 -1 -5 -5 -5 11 11 11 7 7 7 P2 0 -1 -5 11 7 0 - 1 2 3 -4 -6 1 1 1 1 1 1 1 1 P1 0 4 4 -1 14 0 - 0 0 0 - - - 0 -8 -2 10 3 1 2 0 0 0 - - - P2 0 -1 -5 11 7 0 - 0 0 0 - - - 2.1 Gauß-Elimination a pivot b marked 3. Pivotzeile Markieren & Verteilen:wenn Pi Gewinner, dann marked[r]=1, pivot[r]=k, Elemente a[r][k]... a[r][n] und b[r] in Hilfsvektor kopieren und an alle andere Prozessoren senden 4. Berechnung der Eliminationsfaktoren:Pi berechnet für seine lokalen Zeilen i, mit i =1,...n/p, mit marked[i]=0, die Eliminationsfaktoren

  14. 1 2 3 -4 -6 P1 0 4 4 -1 14 0 -8 -2 10 3 P2 1 2 3 -4 -6 1 1 0 -1 -5 11 7 P1 0 0 5 4 15,5 0 - 0 -8 -2 10 3 1 2 1 1 P2 0 0 -4,75 9,75 6,625 0 - 0 - 0 - 0 - 2.1 Gauß-Elimination • 5. Neu Berechnung Matrixelemente:Pi berechnet seine Elemente a[i][j], und b[i], i=1...n/p, j=k+1...n, mit marked[i]=0, neu. marked a b pivot • jeder Prozessor sucht in seinem Array marked nach marked[j]=0 und setzt pivot[j]=n • Lösung des Gleichungssystems durch Rückwärtseinsetzten

  15. Pivotzeile kopieren pivot Eliminationsfaktoren & Elemente von A, b neu berechnen Initialisierung marked Pivotzeile ermitteln Kommunikationsaufwand: 2.1 Gauß-Elimination Berechnungsaufwand:

  16. 2.1 Gauß-Elimination • Vorteile: • Vorhersagbarkeit der Laufzeit, des Speicherbedarfs • Berechnung einer exakten Lösung • Nacheilt: • bei dünnbesetzten Matrizen entsteht fill-in

  17. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

  18. 2.2 zyklische Reduktion • Eine Matrix A mit A= (aji) i,j=1,...,nRnn heißt Bandmatrix mit halber Bandbreite r, falls aij = 0 für |i-j| >r. • Falls r = 1 , so heißt A Tridiagonalmatrix • Gauß-Elimination hier nicht Sinnvoll!

  19. 2.2 zyklische Reduktion • Falls A symmetrisch und positiv definit Lösung mittels • rekursiven Verdoppelns oder zyklischer Reduktion • Ausgangspunkt beider Verfahren: Idee: • xi-1und xi+1aus der i-ten Gleichung, durch einsetzen von geeigneten Vielfachen der Gleichungen i+1 und i-1, zu eliminieren

  20. 2.2 zyklische Reduktion • Nach log n  Schritten, nur noch Elemente der Hauptdiagonalen ungleich Null • Lösung ist einfach abzulesen • Pro Schritt werden O(n) Werte berechnet: • O(nlog n) gegenüber O(n) für Gauß-Elimination aber Berechnung der Werte innerhalb eines Schrittes parallel möglich

  21. 2.2 zyklische Reduktion rekursives Verdoppeln : • Diagonalmatrix in log n Schritten ermitteln zyklische Reduktion : • die Zahl der berechneten Werte in jedem Schritt halbiert • erfordert abschließende Substitutionsphase

  22. 2.2 zyklische Reduktion parallele Implementierung: • Zeilen der Tridiagonalmatrix A blockweise auf p Prozessoren aufgeteilen • Zur Vereinfachung: • n = pq mit q • q = 2Q mit Q  N • Pi speichert Zeilenblock der Größe q, mit den Zeilen (i-1)q+1,...,i*q, für 1  i  p

  23. log p Stufen Q Stufen Q Stufen x1 x2 x3 x4 x5 x6 x7 x8 2.2 zyklische Reduktion parallele Implementierung: 1. zyklische Reduktion:bis nur noch eine Gleichung pro Prozessor vorhanden 2. rekursives Verdoppeln:das nur noch p-dimensionale Gleichungssystem wird nun mittels rekursiven Verdoppeln in gelöst 3. Substitutionsphase:Ermittlung der noch fehlenden Werte des Lösungsvektors x

  24. Aufwand Phase 2: • Aufwand Phase 3: 2.2 zyklische Reduktion • Aufwand Phase 1:

  25. 2.2 zyklische Reduktion • Aufwand Insgesamt:

  26. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

  27. 3.1 klassische iterative Verfahren • Berechnen Folge von Approximationsvektoren{x(k)}k=1,2,...,,die gegen die gesuchte Lösung x*Rn konvergieren • Zerlegung der Matrix A in A=M-N mit M,N Rnn • M nichtsinguläre Matrix • M-1 leicht zu berechnen, z.B. Diagonalmatrix • x* erfüllt Gleichung: Mx* = Nx* +b. • Iterationsvorschrift: Mx(k+1) = Nx(k) + b

  28. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

  29. In Komponentenschreibweise: 3.1.1 Jacobi-Verfahren • Zerlegung von A in A=D-L-R mit (D,L,R Rnn), • D Diagonalmatrix, L untere Dreiecksmatrix, R obere Dreiecksmatrix (jeweils ohne Diagonale) • Iterationsvorschrift: Dx(k+1) =(L+R)x(k) + b

  30. 3.1.1 Jacobi-Verfahren • Abbruchkriterium: relativer Fehler ||.|| Vektornorm, z.B. ||x|| = max i=1,...,n|xi| oder ||x||2=(ni=1|x|2)½ .

  31. 3.1.1 Jacobi-Verfahren parallel Implementierung: • Pimit i=1,...,p speichert die Werte xi(n/p),..., x(i+1)(n/p)-1inklusive dazugehörige Zeilen von A und der Werte von b • jeder Pi führt nun alternierend folgende Aktionen aus, bis Abbruchkriterium erfüllt: 1. Pi sendet seine lokal gehaltenen Elemente des Vektors x(k) an alle anderen Prozessoren (All-to-All Broadcast) 2. Pi berechnetx(k+1) für seine Elemente xj(k+1) mit j=i(n/p),...,(i+1)(n/p)-1 3. Test des Abbruchkriterium

  32. 3.1.1 Jacobi-Verfahren • Aufwand : • In Schritt 1. sendet Pi n/p Werte an p-1 Prozessoren Kommunikationsaufwand Θ((p-1)  (n/p)). • In Schritt 2. n/p Elemente des Vektors x(k+1) berechnen, mit Rechenaufwand von Θ(n (n/p)) • Der Test auf Abbruch in 3. erfordert einen Aufwand von Θ (log p). • Bewertung: • Konvergenz nicht gesichert , niedrige Konvergenzrate • es 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 • kein fill-in, geringer Speicherbedarf bei dünnbesetzten Matrizen

  33. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

  34. 3.1.2 Gauß-Seidel-Verfahren • Zerlegung der Matrix A anlog zum Jacobi-Verfahren, A=D-L-R (D,L,R Rnn) • Iterationsschritt: (D-L)x(k+1)=Rx(k)+b. • in Komponentenschreibweise

  35. 3.1.2 Gauß-Seidel-Verfahren • Einbeziehung der neu Berechneten Approximation • deutlich verbesserte Konvergenzrate • Konvergenzrate: n(e/3)Iterationen um den Fehler um den Faktor 10-e zu senken • aber Datenabhängigkeiten, die eine Parallelisierung des Verfahrens erschweren

  36. Bsp. Temperaturverlauf im Wasserbad xi,j= (xi-1,j + xi+1,j + xi,j-i + x i,j+1) / 4 3.1.2 Gauß-Seidel-Verfahren • In Dünnbesetzte Matrizen weniger Datenabhängigkeiten • Eine Möglichkeit Rot-Schwarz-Anordnung, für Gleichungssysteme, die von einer Gitterstruktur herrühren

  37. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 3.1.2 Gauß-Seidel-Verfahren • Für n = 16 Punkte ergibt sich folgendes Gitter und Matrix A

  38. 1 2 3 4 5 6 7 8 9 10 11 12 1 9 2 10 13 14 15 16 1 2 11 3 12 4 3 4 5 13 6 14 5 6 15 7 16 8 7 8 3.1.2 Gauß-Seidel-Verfahren • Rote-Schwarz-Anordnung : • rot Punkt nur schwarze Nachbarn und umgekehrt • Die Rote von 1,...,nS, numerieren • die Schwarze von nS+1,... ,nR+nS numerieren (nR+nS=n) • neue Matrix  (Permutation von A) • Die Zerlegung von  in  =D-L-R (D,L,R Rnn) mit

  39. 3.1.2 Gauß-Seidel-Verfahren • Iterationsschritt: • in Komponentenschreibweise:

  40. 3.1.2 Gauß-Seidel-Verfahren parallel Implementierung: • maximal p=nr (bzw. p=ns) Prozessoren • Aufteilung anhand des Gitters vorgenommen • wird Pj Gitterpunkt i zugeordnet, dann Pj für Berechnung der Approximationen von xi zuständig • Aufteilung z.B. zeilenblockweise • bei p Prozessoren erhält jeder Prozessor n / p Gitterzeilen (½(n / p) rote und schwarze Punkte)

  41. 3.1.2 Gauß-Seidel-Verfahren • wähle beliebigen Startvektor • Prozessoren iterieren folgende Schritte bis der Test auf Abbruch erfolgreich 1. Pi sendet seine lokal Elemente von xS(k) an alle anderen Prozessoren (All-to-All Broadcast) 2. Pi berechnet die neue Approximation xR(k+1) für seine Elemente von xR(k) 3. Pi sendet seine Elemente von xR(k+1) an alle anderen Prozessoren (All-to-All Broadcast) 4. Pi berechnet die neue Approximation xS(k+1) für seine Elemente von xS(k) 5. Test auf Abbruch

  42. 3.1.2 Gauß-Seidel-Verfahren • Der Rechenaufwand fast identisch zu dem des Jacodi-Verfahrens • In Schritt 2 und 4 jeweils n/2p Elemente der neuen Approximation berechnen, • insgesamt die gleiche Anzahl, wie beim Jacobi-Verfahrens • eine All-to-All Broadcastoperation mehr erforderlich (Schritt 1 und 3) • zusätzlicher Aufwand durch verbessertes Konvergenzverhalten kompensiert Aufwand:

  43. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

  44. 3.2 Methode der konjugierten Gradienten • Voraussetzung: Matrix ARnnpositiv definit und symmetrisch • Lösung x* des Gleichungssystems entspricht dem Minimum der Funktion: • f(x) ist konvex und besitzt eindeutiges Minimum x*

  45. 3.2 Methode der konjugierten Gradienten • Iterationsschritt: • (k) Schrittweite , p(k) Richtungsänderung • Bestimmung von (k) bei bekanntem p(k) : • Notwendige Bedingung, mit r = Ax-b : • Hinreichende Bedingung: immer erfüllt !

  46. 3.2 Methode der konjugierten Gradienten • Wahl des Richtungsvektors p(k) : • die Richtungsvektoren p(k) sind so zu wählen, dass sie konjugierte Basis des Rn bzgl. A bilden: • 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 eine konjugierte Basis des Rn bzgl. A bilden, dann liegt die exakte Lösung nach maximal n Schritten vor.

  47. 3.2 Methode der konjugierten Gradienten Der Algorithmus: • Initialisierung: wähle Startvektor x(0) und setzte p(0) =-r(0) =b-Ax(0) sowie k=0 • Iteration: solange ||r(k)||> berechne 1. q(k) =Ap(k) 2. k = [(r(k))T r(k)] / [ (p(k))T q(k) ] 3. x(k+1) =x(k) + kp(k) 4. r(k+1) =r(k) + kq(k) 5. k= [(r(k+1))T r(k+1)]/ [(r(k))T r(k)] 6. p(k+1) =-r(k+1)+ kp(k) 7. k++ Basisoperationen: Matrix-Vektor-Multiplikation ein inneres Produkt eine Vektor-Addition eine Vektor-Addition ein inneres Produkt eine Vektor-Addition

  48. 3.2 Methode der konjugierten Gradienten parallele Implementierung: • zeilenblockweise Aufteilung von A und blockweise Aufteilung der Vektoren • durch Parallelisierung der verwendeten Basisoperationen • eine Matrix-Vektor-Multiplikation • zwei innere Produkte • drei Vektor-Additionen

  49. 3.2 Methode der konjugierten Gradienten Rechenaufwand:

  50. Gliederung 1. Einleitung 2. Direkte Verfahren • 2.1 Gauß-Elimination • 2.2 zyklische Reduktion 3. Iterative Verfahren • 3.1 klassische iterative Verfahren • 3.1.1 Jacobi-Verfahren • 3.1.2 Gauß-Seidel-Verfahren • 3.2 Methode der konjugierte Gradienten 4. Zusammenfassung

More Related