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:

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:

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:

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