Der eigentliche SGEFA-Algorithmus verwendet eine dreifach geschachtelte Schleife:
DO k = 1, N
find pivot
DO j = k+1, N
call SAXPY(N-k, -A[k,j], A[k+1,k], A[k+1,j],1)
Die naheliegenste Parallelisierung dieses Algorithmus, wie sie zum Beispiel Intel auf dem iPSC/860 verwendete [InScott], sieht so aus:
DO k = 1, N
if (I own the pivot)
broadcast pivot column
else
receive pivot column
DO j over my active columns
call SAXPY
Es wird N mal kommuniziert, der serielle Anteil wächst mit O(N).
Um die Prozessoren gut auszulasten, wurden die Spalten der Matrix nicht fortlaufend auf die Prozessoren verteilt, sonderen nach folgendem Schema (Beispiel: 3 Prozessoren):
Prozessor 1 enthält Elemente 1 und 4, Prozessor 2 Elemente 2 und 5, Prozessor 3 Elemente 3 und 6 (siehe auch Abschnitt 4.3 auf Seite
). Diese Verteilung wird auch Wrapping genannt.
Damit erreichte Intel in [InScott] auf dem iPSC/860 mit 32 Prozessoren 126 MFlop/s. Eine grössere Anzahl Prozessoren brachte keine Performancesteigerung.
Um eine weitere Steigerung zu erreichen, implementierten wir unser SGEFA mit folgenden Änderungen:
Der SGESL Algorithmus hat einen kleineren Rechenaufwand. Er verwendet nur zwei doppelt geschachtelte Schleifen:
DO k = 1,N-1
l=IPVT[k]
t=B[l]
B[l]=B[k]
B[k]=t
call SAXPY(N-k, t, A[k+1,k], 1, B[k+1], 1)
DO k = N,1
B[k]=B[k]/A[k,k]
t=-B[k]
call SAXPY(k-1, t, A[1,k],1,B[1],1)
Eine Parallelisierung dieses Algorithmus lässt sich nur für den SAXPY erreichen:
DO k = 1,N-1
if (I own B[IPVT[k]])
send B[IPVT[k]]
else
receive B[IPVT[k]]
if (I own B[k])
send B[k]
else
receive B[k]
call SAXPY
DO k = N,1
if (I own B[k])
B[k]=B[k]/A[k,k]
send B[k]
else
receive B[k]
call SAXPY
Es sind N Kommunikationen nötig. Der serielle Anteil steigt mit O(N). Beim MUSIC lohnt sich diese Parallelisierung, da ein broadcast to all nicht zeitaufwendiger ist als eine Message nur an einen Adressaten. Auf einem Messagepassingsystem könnte auch diese Variante interessant sein:
Send A to processor 0
if (I am processor 0)
call SGESL
Wir haben hier zwar keinen parallelen Anteil mehr, dafür konnten wir die ineffizienten broadcasts to all umgehen.
Auf einem Einprozessorsystem kann der Aufwand für SGESL gegenüber SGEFA vernachlässigt werden. Bei einem Mehrprozessorsystem ist dies nicht mehr der Fall, da viel kommuniziert werden muss, und der serielle Anteil in der gleichen Grössenordnung wie bei SGEFA liegt. In der Literatur finden sich sehr viele Ansätze über die Parallelisierung der Matrixfaktorisierung, die anschliessende Auflösung hingegen wird selten erwähnt.