Wir haben uns dafür entschieden, den Parallelisierungsschritt auf Stufe LINPACK und LAPACK zu machen. Eine Parallelisierung auf Stufe der BLAS-Routinen würde ihrem Zweck widersprechen: sie wären nicht mehr einfach, universell anwendbar. Der Sinn der BLAS-Routinen ist es aber, möglichst flexibel zu sein. Die BLAS-Routinen werden von verschiedenen LINPACK und LAPACK Routinen aufgerufen, eine Assembleroptimierung einzelner BLAS Routinen bewirkt somit für mehrere Routinen eine Performancesteigerung.
Bei einzelnen BLAS-Routinen haben wir nicht den ganzen Funktionsumfang optimiert. Diese Routinen erkennen aber, ob die optimierte Version aufgerufen werden kann, verwenden sonst aber die allgemeine Routine. Die optimierten Versionen können auch direkt aufgerufen werden, was eine grössere Performance bewirkt. Beim direkten Aufruf werden die vorgegebenen Parameter weggelassen, und der Name der Routine zeigt an, welche Bedingungen erfüllt sein müssen. Für die genaue Syntax siehe Quelltext der C-Versionen.
Das Konzept, dass die BLAS-Routinen nicht angepasst werden hat sich im Falle der ungeblockten Faktorisierung (SGEFA) nicht bewährt. Wir sind deshalb dort von unserem Konzept abgewichen und haben eine entsprechend angepasste SAXPY-Routine programmiert.
Die BLAS-Routinen in C sind allgemein geschrieben und sollten auf jedem System funktionieren. Die Assembler-Versionen nutzen die Fähigkeiten des DSP 96002 aus. Sie gehen in der Regel davon aus, dass das prozessorinterne X-Memory ab Adresse 10 zur Verfügung steht und belegen auch das statische RAM.
Tabelle 8.1 auf Seite 
zeigt im Überblick, welche Routinen portiert worden sind.
Tabelle: Übersicht der portierten BLAS, LINPACK und LAPACK Routinen
ISAMAX gibt den Index des Elementes mit dem grössten Absolutwert zurück.
MINT isamax(MINT N, MFLOAT *SX, MINT INCX)
Die Assemblerversion arbeitet nur mit positivem Inkrement INCX.
SCOPY kopiert den Vektor SX in SY.
void scopy(MINT N, MFLOAT *SX, MINT INCX, MFLOAT *SdY, MINT INCY)
Die Assemblerversion arbeitet nur mit positivem Inkrement INCX und INCY.
SAXPY addiert zum Vektor SY den Vektor SX skaliert um den Faktor SA.
void saxpy(MINT N, MFLOAT SA, MFLOAT *SX, MINT INCX, MFLOAT *SY, MINT INCY)
Dies ist die allgemeine Variante und wurde auch in Assembler programmiert. Es existiert eine spezielle Variante des SAXPY-Algorithmus in Assembler, der die übergeordnete Schleife, wie sie die Routine SGEFA verwendet, enthält. Die spezialisierte Version heisst SAXPY_XYF und arbeitet mit einer anderen Arrayindexierung (siehe 5.1 auf Seite ).
Beim abwechselnden Zugriff auf SX und SY würde mit dem DSP 96002 jedes Mal ein Pagefault entstehen. Damit würden wir in der Schleife von SGEFA 5+5+3=13 Zyklen benötigen. Kopieren wir SX hingegen zuerst in das prozessorinterne X-Memory (3 Zyklen), kommen wir in der Schleife mit 3+3 Zyklen aus. Somit können wir uns von 13 Zyklen auf 9 verbessern, was eine Performanceverbesserung in der Schleife von 6.15 auf 8.89 Mflop/s bedeutet. Bei der spezialisierten Version kann das kopierte SX mehrmals verwendet werden. Somit wächst die maximale Performance gegen 13.33 Mflop/s, falls der prozessorinterne Speicher beliebig gross ist (theoretischer Wert für eine unendliche Matrix).
Die Routinen SAXPY und SAXPY_XYF verwenden den prozessorinternen X-Speicher ab Adresse 10. Ist der Vektor länger als 502 Byte, kommt der Rest in den statischen Speicher zu liegen. Für diesen Teil liegt nun der Zugriff bei 2+2+5=9 Zyklen. Das Umkopieren der Daten benötigt 2+5 Zyklen, womit wir bei 16 Zyklen liegen. Hier lohnt sich also das Umkopieren nicht. Bei der spezialisierten Routine SAXPY_XYF kann wiederum für grosse Matrizen die Kopierzeit vernachlässigt werden, womit die maximale Performance gegen 8.89 Mflops/s geht.
SSCAL skaliert einen Vektor mit einer Konstante.
void sscal(MINT N, MFLOAT SA, MFLOAT *SX, MINT INCX)
Die Assemblerversion arbeitet nur mit positivem Inkrement INCX.
SSWAP vertauscht die Vektoren SX und SY.
void sswap(MINT N, MFLOAT *SX, MINT INCX, MFLOAT *SY, MINT INCY)
SGER führt die Rank 1 Operation durch: 
void sger(MINT M, MINT N, MFLOAT ALPHA, MFLOAT *SX, MINT INCX, MFLOAT *SY, MINT INCY, MFLOAT *A, MINT LDA)
Die Assemblerversion arbeitet nur mit positivem Inkrement INCX und INCY.
SGEMV führt eine Matrix-Vektor Multiplikation durch: 
void sgemv(MINT TRANSA, MINT M, MINT N, MFLOAT ALPHA, MFLOAT *A, MINT LDA, MFLOAT *X, MINT INCX, MFLOAT BETA, MFLOAT *Y, MINT INCY)
Die Assemblerversion arbeitet nur mit positivem Inkrement INCX und INCY sowie nicht transponiertem A.
SGEMM führt eine Matrix-Matrix Multiplikation durch: 
void sgemm(MINT TRANSA, MINT TRANSB, MINT M, MINT N, MINT K, MFLOAT ALPHA, MFLOAT *A, MINT LDA, MFLOAT *B, MINT LDB, MFLOAT BETA, MFLOAT *C, MINT LDC)
Wenn A und B nicht transponiert sind, kann die Assemblerversion aufgerufen werden. Die Assemblerversion benützt den prozessorinternen X-Speicher sowie den Beginn des statischen RAMs.
STRSM löst eine der Matrix-Gleichungen: 
oder 
.
void strsm(MINT SIDE, MINT UPLO, MINT TRANSA, MINT DIAG, MINT M, MINT N,
MFLOAT ALPHA, MFLOAT *A, MINT LDA, MFLOAT *B, MINT LDB)
Für Left, Upper, Notranspose, Diagnonunit, Alpha = 1.0 und Left, Lower, Notranspose, Diagunit, Alpha = 1.0 können die entsprechenden Assemblerversionen aufgerufen werden. Die Assemblerversionen benützen den prozessorinternen X-Speicher sowie den Beginn des statischen RAMs.
Diese Routinen verwenden eine andere Arrayindexierung als in Fortran üblich. Bei diesen Routinen ist das nächste Element im Speicher dasjenige, das auf der gleichen Zeile in der nächsten Spalte liegt.
allocate_matrix reserviert den benötigten Speicher für die Faktorisierung und Lösung.
MINT allocate_matrix(MINT N)
matrix_to_music schreibt die Matrix, deren Speicherplatz zuvor mit allocate_matrix reserviert wurde, zum MUSIC.
MINT matrix_to_music(void)
Das Schreiben der Daten zum MUSIC muss wie folgt aussehen:
MCALL ( matrix_to_music() );
Wr_to_music (A, N*DimZ*NumDSPs*sizeof(MFLOAT));
Wr_to_music (B, DimZ*NumDSPs*sizeof(struct Pack));
matrix_from_music liest die Matrix wieder zum Host zurück.
MINT matrix_from_music(void)
Das Lesen der Daten vom MUSIC muss wie folgt aussehen:
MCALL ( matrix_to_music() );
Rd_from_music(A,N*DimZ*NumDSPs*sizeof(MFLOAT));
time_from_music liest die Zeit, die für die Faktorisierung und die Auflösung benötigt wurde, zum Host zurück.
MINT time_from_music(void)
Das Lesen der Zeit muss wie folgt aussehen:
MCALL ( matrix_to_music() );
for (j=0;j<DimZ;j++)
Rd_from_music(&timerec, 3*sizeof(MFLOAT));
wobei
struct time 
MFLOAT fac;
MFLOAT sol;
MFLOAT special;
timerec;
mit timerec.fac = Zeit Faktorisierung (Sekunden)
und timerec.sol = Zeit Faktorisierung (Sekunden)
result_from_music liest den Resultatvektor zum Host zurück.
MINT result_from_music(void)
Das Resultat wird wie folgt gelesen:
MCALL(result_from_music());
Rd_from_music(B, N*sizeof(MFLOAT));
ipvt_from_music liest die Pivotindices zum Host zurück.
MINT ipvt_from_music(void)
Die Pivotindices werden wie folgt gelesen:
MCALL(ipvt_from_music());
Rd_from_music(IPVT, N*sizeof(MFLOAT));
free_matrix gibt allen Speicherplatz wieder frei.
sgefa faktorisiert die Matrix.
MINT sgefa(void)
Vorgängig müssen folgende Routinen aufgerufen worden sein:
sgesl löst die faktorisierte Matrix auf.
MINT sgesl(void)
Vorgängig müssen folgende Routinen aufgerufen worden sein:
SGETF2 berechnet die LU Faktorisierung einer allgemeinen 
Matrix mit partial pivoting und Zeilentauschen.
MINT sgetf2(MINT M, MINT N, MFLOAT *A, MINT LDA, MINT *IPIV
Diese Routine ist in ein eigenes Modul ausgelagert, damit sie optimal im Programmspeicher plaziert werden kann.
SLASWP führt eine Anzahl Zeilenvertauschungen durch.
void slaswp(MINT N, MFLOAT *A, MINT LDA, MINT K1, MINT K2, MINT *IPIV,
MINT INCX)
SGETRF berechnet die LU Faktorisierung einer allgemeinen Matrix blockweise.
MINT sgetrf(MINT M, MINT N, MFLOAT *A, MINT LDA, MINT *IPIV)
SGETRS löst ein der System 
oder 
, wobei A durch SGETRF berechnet worden sein muss.
MINNT sgetrs(MINT TRANS, MINT N, MINT NRHS, MFLOAT *A, MINT LDA,
MINT *IPIV, MFLOAT *B, MINT LDB)
Diese Routinen verwenden die in Fortran übliche Arrayindexierung: das nächste im Speicher liegende Element ist dasjenige, das in der gleichen Spalte eine Zeile tiefer liegt.
Send_ilaenv_para sendet die Blockgrösse vom Host an den MUSIC. Diese Funktion wird im Hostprogramm aufgerufen.
void Send_ilaenv_para(MINT N)
Send_ilaenv_para überprüft, ob die übergebene Blockgrösse möglich ist. Ist sie das, wird sie dem MUSIC mitgeteilt, andernfalls wird ein Default-Wert gesendet.
allocate_matrix reserviert den benötigten Speicher für die Faktorisierung und Lösung.
MINT allocate_matrix(MINT N)
matrix_to_music schreibt die Matrix, deren Speicherplatz zuvor mit allocate_matrix reserviert wurde, zum MUSIC.
MINT matrix_to_music(void)
Das Schreiben der Daten zum MUSIC muss wie folgt aussehen:
MCALL ( matrix_to_music() );
for (j=0;j<DimZ;j++)
Wr_to_music (&A[j*N*BlockSize*NumDSPs],
N*BlockSize*NumDSPs*sizeof(MFLOAT));
Wr_to_music(B, N*sizeof(MFLOAT));
mit DimZ = ((N-1)/Blocksize)/NumDSPs + 1
matrix_from_music liest die Matrix wieder zum Host zurück.
MINT matrix_from_music(void)
Das Lesen der Daten vom MUSIC muss wie folgt aussehen:
MCALL ( matrix_to_music() );
for (j=0;j<DimZ;j++)
Rd_from_music (&A[j*N*BlockSize*NumDSPs],
N*BlockSize*NumDSPs*sizeof(MFLOAT));
mit DimZ = ((N-1)/Blocksize)/NumDSPs + 1
time_from_music liest die Zeit, die für die Faktorisierung und die Auflösung benötigt wurde, zum Host zurück.
MINT time_from_music(void)
Das Lesen der Zeit muss wie folgt aussehen:
MCALL ( matrix_to_music() );
for (j=0;j<DimZ;j++)
Rd_from_music(&timerec, 3*sizeof(MFLOAT));
wobei
struct time
MFLOAT fac;
MFLOAT sol;
MFLOAT special;
timerec;
mit timerec.fac = Zeit Faktorisierung (Sekunden)
und timerec.sol = Zeit Faktorisierung (Sekunden)
result_from_music liest den Resultatvektor zum Host zurück.
MINT result_from_music(void)
Das Resultat wird wie folgt gelesen:
MCALL(result_from_music());
Rd_from_music(B, N*sizeof(MFLOAT));
ipvt_from_music liest die Pivotindices zum Host zurück.
MINT ipvt_from_music(void)
Die Pivotindices werden wie folgt gelesen:
MCALL(ipvt_from_music());
Rd_from_music(IPVT, N*sizeof(MFLOAT));
free_matrix gibt allen Speicherplatz wieder frei.
MINT free_matrix(void)
MUSIC_sgetrf faktorisiert die Matrix.
MINT MUSIC_sgetrf(void)
Vorgängig müssen folgende Routinen aufgerufen worden sein:
MUSIC_sgetrs löst die faktorisierte Matrix auf.
MINT MUSIC_sgetrs(void)
Vorgängig müssen folgende Routinen aufgerufen worden sein:
Das Demoprogramm löst die Gleichung Ax=b. A und b werden ab Disk gelesen, x anschliessend wieder auf Disk geschrieben. Das Programm rechnet den relativen Fehler nach [LinUser] und gibt ihn sowie die erreichte Performance aus. Die Anzahl Floatingpointoperationen werden dabei wie in [Perf] verlangt unabhängig vom verwendeten Algorithmus berechnet:
Das Programm 1000s kann über anonymous ftp in netlib/benchmark als 1000s bezogen werden. Es generiert eine Testmatrix, ruft die Prozeduren zum Lösen auf. Anschliessend wird wieder dieselbe Testmatrix generiert und geprüft, ob das Resultat stimmt. Nach [Perf] wird das Residual wie folgt berechnet:
Nach [LinUser] sei 
das Maximum der Beträge der einzelnen Spaltenvektoren von A, 
der Betrag des Vektors, das heisst die Wurzel aus der Summe der Quadrate der einzelnen Elemente. Im Programm werden diese beiden Ausdrücke einfach durch das maximale Element ersetzt, zusätzlich wird durch die Rechengenauigkeit dividiert, was natürlich ein ganz anderes Resultat ergibt. Mit unserem Testprogramm das sich nach [LinUser] richtet, erhalten wir einen relativen Fehler von ca. 
, mit dem Programm 1000s erhalten wir etwa 10.
Interessant ist das File MUSIC.c. Dies ist das C-Programm, das die Aufrufe des Fortran-Programms MUSIC gerecht umändert. Beim Block-Algorithmus (SGETRF) ist dies nur die Anpassung an die andere Parameterübergabe, beim ungeblockten Algorithmus muss auch die Arrayindexierung angepasst werden. Siehe auch Abschnitt 8.8.
Die Portierung von Fortran nach C ist im Prinzip einfach, wenn man einige Dinge beachtet:
Die Benchmark-Programme von J.Dongarra [Perf] verwenden in ihrer Hauptroutine Fortran. Dieser Fortranquellcode darf nicht abgeändert werden. Wir müssen also die MUSIC-Routinen aus Fortran aufrufen können. Wir haben deshalb eine kleine Routine geschrieben, die von Fortran aufgerufen wird und dann die nötigen MUSIC-Aufrufe vornimmt. Das Programm läuft auf dem Host (Sun). Dabei gilt es folgendes zu beachten: