Kommunikation und Speicherzugriff ist in vielen Mehrprozessorsystemen langsam im Vergleich zur Fliesskommarechenleistung. Durch die Verwendung von Matrix-Matrixoperationen kann einerseits die Kommunikation minimiert und anderseits der Speicherzugriff optimiert werden. Welche Blockgrösse die beste Performance liefert, hängt stark von der Kommunikationsarchitektur, dem Speicheraufbau und der Anzahl Prozessoren ab.
Der SGETRS-Algorithmus zum Faktorisieren einer Matrix sieht so aus:
DO j = 1,N,BlockSize
jb = MIN(N-j+1, BlockSize)
call SGETF2(N-j+1, jb, A[j,j], N, IPVT[j], iinfo)
adjust info and IPVT
call SLASWP for columns 1:j-1
if (j+jb<N)
call SLASWP for columns j+jb:N
call STRSM for the triangular solve
call SGEMM to update the trailing submatrix
Als Parallelisierungsmöglichkeit verwendeten wir denselben Algorithmus wie Intel in [InBlock]. Das Wrapping erfolgt blockweise. Die Verteilung der Matrix erfolgt spaltenweise. Dies hat den Nachteil, dass für das Auflösen der Matrix (SGETRS) keine Parallelisierung mehr möglich ist (siehe auch Abschnitt 4.3 auf Seite 
). Das Pipelining mit dem workin und workout array ermöglicht es, den seriellen Anteil möglichst klein zu halten und die Kommunikation möglichst parallel zur Rechnung ablaufen zu lassen.
Bei diesem Algorithmus wird mit wachsender Blockgrösse die Anzahl Kommunikationen kleiner, bei gleichbleibender total zu kommunizierender Datenmenge. Dies ist vor allem auf Systemen grosser Latenzzeit interessant. Auf der anderen Seite wird die optimale Blockgrösse durch die Speicherarchitektur begrenzt. Die Blöcke sollten z.B im Cache gehalten werden können oder nicht über mehrere Pages verteilt sein. Bei unserer Implementation auf dem MUSIC spielt die Anzahl Kommunikationen nicht so eine grosse Rolle, da das MUSIC-System eine kleine Latenzzeit aufweist, die erst noch unabhängig von der Anzahl Adressaten ist. Auf dem MUSIC fanden wir die optimale Blockgrösse um 4 mit Radius 2 herum. Auf dem Intel iPSC/860 ist die Blockgrösse um 8 mit Radius 4 verteilt [InBlock].
Der verwendete Algorithmus sieht nun so aus:
if (I own block 0)
call SGETF2 to factor first panel
ship factored panel and pivots in workout array
DO i = 1,N, BlockSize
if (I own current block)
copy workout array into workin
else
receive panel and pivots into workin
apply pivots
if (there is work left)
if (I am the next to factorize)
call STRSM on first remaining panel
call SGEMM on first remaining panel
call SGETF2 on that panel
and ship it in workout to other processors
call STRSM on workin array
call SGEMM on workin array
Der Lösungsalgorithmus SGETRS baut auf zwei BLAS-Routinen auf:
call SLASWP apply row interchanges
call STRSM to solve L*Y=B
call STRSM to solve U*X=Y
Ist die Matrix nun spaltenweise verteilt, wäre der Aufwand sie umzuverteilen grösser (DimZ-Kommunikationen + Speicherbewegungen) als sie auf einem Prozessor zu lösen. Da aber die ganze Matrix nicht auf einem Prozessor Platz hat, haben wir den untenstehenden Algorithmus gewählt:
call SLASWP apply row interchanges
DO j = 1,N,BlockSize
jb=MIN(N-j,BlockSize)
if (I own current panel)
call STRSM to solve L*Y=B
if (there is work left)
call SGEMV
send x to other processors
else
receive x
DO j = N,1,BlockSize
if (I own current panel)
call STRSM to solve U*X=Y
if (there is work left)
call SGEMV
send x to other processors
else
receive x
Die Anzahl Kommunikationen nehmen mit der Blockgrösse ab. Dies hat nun einen grossen Einfluss, da die Kommunikation nicht mehr parallel zur Rechnung erfolgt. Der Algorithmus ist auch wesentlich komplizierter geworden und befriedigt nicht mehr. Abhilfe haben wir jedoch nicht gefunden.