Das Verfahren der konjugierten Gradienten (CG-Verfahren) zur Lösung des linearen Gleichungssystems konsumiert bei der Matrix-Vektor-Multiplikation, die innerhalb der Iterationsschleife liegt, einen Großteil seiner Rechenzeit. Es soll im weiteren Verlauf gezeigt werden, wie trotz der aufwendigeren Matrixstruktur diese Multiplikation effizient implementiert werden kann.
Zunächst soll die Struktur eines MCSR-(Modified Compressed Sparse Row)-Matrixformates an einem Beispiel erklärt werden.
Als Beispiel wird die symmetrische Matrix mit der Struktur
für die weiteren Betrachtungen auf eine untere Dreiecksmatrix und die Diagonaleinträge reduziert.
Die Diagonaleinträge nehmen in der Matrix eine Sonderstellung ein, da bei positiv definiten Systemen die Diagonale immer voll besetzt ist und die Matrix oft so skaliert wird, daß die Diagonaleinträge zu eins werden. Da bei Vorwärtseliminationen und Rücksubstitutionen nach Abarbeitung einer Zeile durch den Diagonaleintrag dividiert wird ((4.13),(4.15)), erweist es sich als günstig, die Diagonaleinträge in einem gesonderten Feld
abzuspeichern. Die restlichen Einträge der unteren Dreiecksmatrix, die ungleich Null sind, werden in einem eindimensionalen Feld
abgespeichert. Ein paralleles Indexfeld
gibt an, zu welcher Spalte der jeweilige Eintrag in der Matrix zuzuordnen ist. Ein zusätzliches Feld
legt den Beginn einer neuen Zeile mit einer Referenz in fest, wobei der
letzte Eintrag von
den Beginn der virtuellen
Zeile der Matrix
beschreibt, um das letzte Zeilenende festzusetzen. Die Gesamtlänge der
Felder
und
ist also
.
Der Zugriff auf einen gewissen Matrixeintrag ist mit einer Suche im Indexfeld
verbunden. Möchte man zum Beispiel den Eintrag
ermitteln,
so liefert
für die fünfte Zeile den Spaltenstart 3 und das Spaltenende 5.
Nun sucht man das Indexfeld
von der Stelle 3-5 nach dem Eintrag 3 ab und
findet diesen an der Stelle 4. Damit ist der gewünschte Koeffizient im
Feld
an der Stelle 4 mit dem Wert 5.3 gefunden.
Die Einträge in und
können innerhalb einer
Matrixzeile unsortiert vorliegen. Sortiert man aber nach den Spaltenindizes
, so kann man mit einer binären Suche den Aufwand reduzieren.
Die Anzahl der Einträge der Zeile
soll mit
bezeichnet werden, und
die Anzahl der Einträge in die untere Dreiecksmatrix
.
Da bei vielen Matrixstrukturen unter Verwendung quadratischer Ansatzfunktionen
die gemittelte Anzahl der Einträge pro Zeile von
bei etwa
liegt, findet man mit maximal vier Teilungen das Auslangen.
Veranschlagt man für das Gleitkommafeld acht Bytes pro Eintrag und
für das Indexfeld
vier Bytes, so ergibt sich unter der Annahme einer
vollbesetzten Matrix ein zusätzlicher Speicheraufwand von etwas mehr als
gegenüber einer unkomprimierten Matrix. Da die gemittelte Anzahl
der Einträge pro Zeile nicht mit der Größe des Ranges korreliert, wird
die Effizienz des MCSR-Verfahrens um so besser, je größer der Rang des
Gleichungssystems ist.
Die Matrix-Vektor-Multiplikation ist die einzige Operation, die bei der
CG-Methode direkt auf die Matrix angewendet wird. Da für die Multiplikation
die Gesamtmatrix benötigt wird, jedoch nur die untere triangulare Matrix und Diagonale
vorhanden ist, soll die
Multiplikation in die Teilaufgaben
zerlegt werden.
Den Multiplikationsanteil der Diagonale und der unteren Dreiecksmatrix, ausgehend von dem MCSR-Matrix-Format, zeigt der folgende Pseudo-Code:
Den oberen Matrixanteil berechnet man ebenso einfach durch das
Programm:
Bemerkenswert ist, daß keinerlei Suche von Matrixeinträgen notwendig ist und Multiplikationen mit Einträgen, die Null sind, vollständig unterdrückt werden.
Das bedeutet, daß trotz eines Zeilenkompressionsverfahrens der Anteil
ohne Mehraufwand
gegenüber dem Anteil
berechnet werden kann.
Die Programmteile (4.3) und (4.4) lassen sich natürlich zu
zusammenfassen.
Der Aufwand der Vektor-Matrix-Multiplikation liegt, durch das MCSR-Format bedingt,
statt bei multiply-adds bei
multiply-adds.
Für den in Abschnitt 4.5.4.1 gezeigten ICCG-Gleichungslöser ist pro Iteration eine Vorwärtselimination und Rücksubstitution, die ebenso auf dem MCSR-Format arbeitet, durchzuführen. Die Implementierung der Vorwärtselimination und Rücksubstitution hat große Ähnlichkeit mit den gezeigten Teilmultiplikationen (4.3) und (4.4), und ist daher genauso effizient machbar.