1. Bestimme obere Dreiecksmatrix U mit U TU = A .zu 1.)
2. Bestimme y mit U T · y = b .
3. Bestimme x mit Ux = y .
Dann gilt Ax = b , denn aus Ux = y folgt U TUx = U Ty mit U TU = A und U Ty = b .
Sei U bis zur Zeile i - 1 bereits bestimmt. Aus
aii = (uki)2aij = (uikT · ukj)
Zeilenweises Bestimmen der Matrix U | |
zunächst der dunkle, dann der helle Teil |
FOR i:= 0 TO n-1 DO (* i-te Zeile *) tmp := a[i,i]; FOR k := 0 TO i-1 DO tmp := tmp - u[k,i]*u[k,i] END; u[i,i] := sqrt(tmp); (* Diagonalelement *) FOR j := i+1 TO n-1 DO tmp := a[i,j]; FOR k := 0 TO i-1 DO tmp := tmp - u[k,i]*u[k,j]; END u[i,j] := tmp/u[i,i] END END
zu 2.)
bi = uijT · yj
yi : = (bi - uji · yj)/uii
FOR i := 0 TO n-1 DO tmp := b[i]; FOR j := 0 TO i-1 DO tmp := tmp - u[j,i]*y[j] END; y[i] := tmp/u[i,i] END
zu 3.)
yi = uij · xj
xi : = (yi - uij · xj)/uii
FOR i := n-1 DOWNTO 0 DO tmp := y[i]; FOR j := i + 1 TO n-1 DO tmp := tmp - u[i,j]*x[j] END; x[i] := tmp/u[i,i] END
Zur parallelen Cholesky-Zerlegung wird ein Ring von n Prozessoren verwendet. Zu Beginn speichert Prozessor j Spalte j von Matrix A . Während der Rechnung ermittelt Prozessor j die Spalte j von Matrix U . Dabei können die Matrixelemente von A mit denen von U überschrieben werden. Für die parallele Backward-Substitution ist es erforderlich, daß Prozessor i über die i -te Zeile von U verfügt. Dies kann dadurch erreicht werden, daß während der parallelen Zerlegungsphase die entsprechenden Matrixelemente beim Durchreichen einbehalten werden (nicht im Algorithmus berücksichtigt!). Die parallele Forward-Substitution kann vereinfacht werden, indem b als zusätzliche Spalte n von A schon in der Zerlegung behandelt wird (nicht im Algorithmus berücksichtigt).