The Art of Designing Matrix Classes
Natürlich ist die Teilnahme an den Übungen freiwillig. Alternativ zum Votieren im E44 kann auch das submit-Kommando benutzt werden. Hinweise dazu auf der Seite 5.
Wir werfen zunächst einen kritischen Blick auf die low-level BLAS Funktionen im Namensraum hpc::ulmblas. Hier werden Matrizen ähnlich wie in C behandelt:
-
Eine Matrix wird durch die Dimensionen, einen Zeiger auf das erste Element und den Inkrementen beschrieben:
typename <typename Index, typename Alpha, typename TX> void gescal(Index m, Index n, const Alpha &alpha, TX *X, Index incRowX, Index incColX) { /* ...*/ }
Funktionen haben deshlab eine recht große Anzahl an Parametern.
-
Der Zugriff auf Elemente erfolgt durch Pointer Arithmetik:
/* ... */ X[i*incRowX+j*incColX] *= alpha; /* ... */
Dies verlangt ein Verständnis dafür, wie Matrizen im Speicher liegen. Der Programmierer muss sich zudem sicher sein, dass nur auf zulässige Elemente mit \(0\leq i < m\) und \(0\leq j < n\) zugegriffen wird. Copy-Paste-Fehler, die zum Beispiel zu Fehlern wie incRowX statt incColX führen, können nur durch Sorgfalt vermieden oder aufgespürt werden.
-
Die Speicherverwaltung liegt im Verantwortungsbereich des Programmierer
void foo() { double *A = new double[5*3]; /* ... */ hpc::ulmblas::gescal(5, 3, 1,5, A, 1, 5); /* ... */ delete [] A; // nicht vergessen! }
Diese Zusammenfassung soll allerdings keine generelle Kritik darstellen. Ein Vorteil ist, dass nur eine einfache Teilmenge der C++ Programmiersprache verwendet wird. Im Prinzip sind Programmierkenntnisse in C ausreichend.
Für die Implementierung von numerischen Verfahren kann es aber vorteilhaft sein, wenn diese technischen Aspekte abstrahiert werden. Dazu wurden im Namensraum hpc::matvec Matrix-Klassen eingeführt. Die wesentlichen Eigenschaften dieser Klassen fassen wir zunächst zusammen. Anschliessend überlegen wir uns, ob es möglich wäre, dies einfacher zu realisieren.
Ein Blick auf das Vorlesungsprojekt
Im Vorlesungsprojekt der letzten Session wurden für eine rechteckige Matrix (general matrix with full storage) drei verschiedene Klassen implementiert:
-
Eine Instanz/Objekt der Klasse GeMatrix verwaltet einen zugehörigen Speicherbereich. Im Konstruktor wird dieser allokiert und im Destruktor wieder freigegeben.
GeMatrix<double> A(5, 6); // Allokiert Speicher fuer eine 5x6 Matrix
-
Eine Instanz der Klasse GeMatrixView verwendet lediglich einen bereits allokierten Speicherbereich. Dieser kann zum Beispiel von einer GeMatrix allokiert worden sein:
// Matrix A allokiert Speicher fuer eine 5x6 Matrix GeMatrix<double> A(5, 6); // Matrix B referenziert einen Teil der Matrix A: 2 Zeilen und 3 Spalten // beginnend bei A.data (der Adresse von A(0,0)). GeMatrixView<double> B(2, 3, A.data, A.incRow, A.incCol); // Analog haette man diese View wie folgt erzeugen koennen: GeMatrixView<double> C(2, 3, &A(0,0), A.incRow, A.incCol);
Ebenso kann aber auch ein beliebiger Speicherbereich verwendet werden. Es muss lediglich sichergestellt sein, dass dieser große genug ist:
double A[5*6]; // Array auf Stack oder Datensegment /* ... */ double B = new double[5*6]; // Array auf Heap // Views die den Speicher der obige Arrays verwenden GeMatrixView<double> C(2, 3, A, 2, 3); GeMatrixView<double> D(2, 3, B, 2, 3);
Der häufigste Anwendungsfall für Views wird aber das Referenzieren von Matrix-Blöcken darstellen. Um dies zu vereinfachen wurde in allen Matrix-Klassen der Funktionsoperator überladen. Zum Beispiel in GeMatrix:
template <typename T, typename I=std::size_t> struct GeMatrix { /* ... */ typedef GeMatrixView<T,Index> View; /* ... */ View operator()(Index i, Index j, Index m, Index n) { assert(i+m<=numRows); assert(j+n<=numCols); return View(m, n, &(operator()(i,j)), incRow, incCol); } /* ... */ };
Die Implementierung der GeMatrixView Klasse unterscheidet sich von der GeMatrix Klasse im wesentlichen nur bei den Konstruktoren und Destruktoren:
-
In den Konstruktoren wird kein Speicher allokiert. Statt dessen wird der Daten-Zeiger auf bereits allokierten Speicher gesetzt.
-
Der Destruktor tut nichts: Es wurde kein Speicher allokiert und deshalb muss auch keiner freigegeben werden.
-
-
Mit GeMatrixConstView wurde eine dritte Klasse für Matrizen implementiert. Diese gleicht im wesentlichen wiederum der Klasse GeMatrix View. Allerdings wurden alle Methoden entfernt, die es ermöglichen wurde die Matrix zu verändern. Zusätzlich wurde der Daten-Zeiger als const definiert.
Eventuelle Nachteile dieser Realisierung
Für einen numerischen Algorithmus soll es keine Rolle spielen, ob eine Matrix ihren eigenen Speicher verwaltet oder nicht. Um an dieser Stelle nicht auf numerische Fragestellungen eingehen zu müssen, betrachten wir aber nur eine Funktion print zur Ausgabe einer Matrix. Wir benötigen bei unserem Ansatz drei Versionen, die sich nur in der Signatur unterscheiden:
-
Ausgabe einer GeMatrix:
template <typename Index, typename T> void print(const GeMatrix<T, Index> &A, const char *name = "") { if (*name) { printf("%s = \n", name); } for (Index i=0; i<A.numRows; ++i) { for (Index j=0; j<A.numCols; ++j) { print_value(A(i,j)); } printf("\n"); } printf("\n"); }
-
Ausgabe einer GeMatrixView:
template <typename Index, typename T> void print(const GeMatrixView<T, Index> &A, const char *name = "") { if (*name) { printf("%s = \n", name); } for (Index i=0; i<A.numRows; ++i) { for (Index j=0; j<A.numCols; ++j) { print_value(A(i,j)); } printf("\n"); } printf("\n"); }
-
Ausgabe einer GeMatrixView:
template <typename Index, typename T> void print(const GeMatrixConstView<T, Index> &A, const char *name = "") { if (*name) { printf("%s = \n", name); } for (Index i=0; i<A.numRows; ++i) { for (Index j=0; j<A.numCols; ++j) { print_value(A(i,j)); } printf("\n"); } printf("\n"); }
Für Funktionen mit weiteren Matrix-Argumenten explodiert die Anzahl der Varianten. Störend und fehleranfällig ist zudem die Implementierung im Stil von Copy-Paste. Und noch schlimmer, ändert man anschliessend eine Variante, so muss man analog alle anderen Varianten ändern.
Verwendung von Template-Funktionen
Diese stupide Arbeit können wir durch Templates dem Compiler überlassen:
template <typename Matrix> void print(const Matrix &A, const char *name = "") { typedef typename Matrix::Index Index; if (*name) { printf("%s = \n", name); } for (Index i=0; i<A.numRows; ++i) { for (Index j=0; j<A.numCols; ++j) { print_value(A(i,j)); } printf("\n"); } printf("\n"); }
Dies funktioniert wunderbar, solange man keine print Funktion für andere Datentypen im gleichen Namensraum definieren möchte. Stellen wir uns vor, dass wir irgendwann auch für Vektoren entsprechende Klassen implementieren. Dann wäre es naheliegend eine zugehörige print Funktion zu implementieren:
template <typename Vector> void print(const Vector &x, const char *name = "") { typedef typename Vector::Index Index; if (*name) { printf("%s = \n", name); } for (Index i=0; i<x.length; ++i) { print_value(x(i)); } printf("\n"); }
Bei einem Aufruf wie in
void foo() { GeMatrix<double> A(4,5); /* ... */ print(A); }
wäre für den Compiler nicht entscheidbar, welcher Funktionsaufruf erfolgen soll. Ein Lösung wäre es den Typ der Objekte in den Funktionsnamen aufzunehmen. Dies würde zu printGeMatrix und printVector führen. Bei Funktionen mit mehreren Matrix/Vector Parametern würde dies im Allgemeinen zu sehr langen Funktionsnamen führen. Auch diese stupide Erzeugung von unterschiedlichen Funktionsnamen können wir dem Compiler überlassen. Um festzulegen welche Funktion aufgerufen werden soll verwenden wir das SFINAE Konzept mit entsprechenden Traits:
template <typename MA> typename std::enable_if<IsGeMatrix<MA>::value, void>::type print(const MA &A, const char *name = "") { typedef typename MA::Index Index; if (*name) { printf("%s = \n", name); } for (Index i=0; i<A.numRows; ++i) { for (Index j=0; j<A.numCols; ++j) { print_value(A(i,j)); } printf("\n"); } printf("\n"); }
Für jede Matrix/Vector Klasse muss dann aber ein entsprechender Trait wie hier IsGeMatrix implementiert werden:
template <typename Any> struct IsGeMatrix { static constexpr bool value = false; }; /* ... */
Aufgaben
Implementiert zwei Template-Klassen Foo und Dummy und zugehörige Traits, so dass dieses Beispiel
#include <cstdio> #include <type_traits> /* ... */ template <typename T> typename std::enable_if<IsFoo<T>::value, void>::type print(const T &) { std::printf("foo\n"); } template <typename T> typename std::enable_if<IsDummy<T>::value, void>::type print(const T &) { std::printf("dummy\n"); } int main() { Foo<double> foo; Dummy<float> dummy; print(foo); print(dummy); }
Diese Ausgabe erzeugt:
foo dummy