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.