C++-Programmierung/ Eine Matrix-Bibliothek – mitrax/ Operationen

In diesem Kapitel wollen wir nun endlich die Rechenoperationen einführen, die für eine mitrax-Matrix standardmäßig definiert sein sollen. Dies sind die Addition und Subtraktion zweier Matrizen, die Skalarmultiplikation und die Matrixmultiplikation, sowie für quadratische Matrizen die Berechnung der Determinante und die Bildung der inversen Matrix, sofern die Matrix regulär ist. Außerdem wollen wir eine Negation einführen und natürlich sollten zwei Matrizen miteinander verglichen werden können.

Sofern nicht anders angegeben, werden die Operationen in diesem Kapitel in der Datei simple_operator.hpp definiert. „Simple“ deshalb, weil wir im nächsten Kapitel noch Operatoren auf Basis einer verzögerten Berechnung einführen wollen.

Negation

Das negieren einer Matrix erfolgt Elementweise. Die Implementierung kann somit hervorragend über das Iteratorinterface der Matrix erfolgen. Wie verwenden den uniären Minusoperator um eine Negation durchzuführen. Dabei wird ein neues Objekt erzeugt. Wenn Sie eine Zahl negieren, ist das Ergebnis schließlich auch eine neue Zahl. Die Verwendung des uniären Minusoperators schürt die Erwartung, dass auch ein uniärer Plusoperator existieren muss. Auch dies ist äquivalent zu gewöhnlichen Zahlen. Sie können ein Plus davor schreiben, die Zahl verändert sich aber nicht. Dennoch erzeugt auch dies ein neues Objekt. Dieser Plusoperator dient einzig als übliche Ergänzung zum Negationsoperator. Praktischen nutzen hat er im Grunde keinen, außer dass der Gegensatz zu einer Negation, an bestimmten Stellen, deutlicher ausgedrückt wird, wenn das Plus davor explizit angegeben ist.

template < typename T, typename Container >
inline matrix< T, Container > const operator-(matrix< T, Container > op){
    for(typename matrix< T, Container >::iterator i = op.begin(); i != op.end(); ++i){
        *i = -*i;
    }
    return op;
}

template < typename T, typename Container >
inline matrix< T, Container > const operator+(matrix< T, Container > const& op){
    return op;
}

Addition und Subtraktion

Die Addition und Subtraktion erfolgt für Matrizen ebenfalls elementweise. Daher müssen zwei Matrizen, die in dieser Weise miteinander Verknüpft werden sollen, die gleiche Dimension haben. Folglich müssen wir ein error::dimension_unequal-Objekt werfen, wenn diese Bedingung nicht erfüllt ist. Beide Operationen werden im entsprechenden verketteten Zuweisungsoperator implementiert. Die gewöhnlichen Operatoren rufen ihre mit der Zuweisung verketten Äquivalente auf, indem sie zunächst den ersten Parameter kopieren und anschließend für das neue Objekt, die verkettete Operation ausführen. Den Minusoperator könnten wir prinzipiell auch durch eine Addition mit der negierten Matrix implementieren, dies würde jedoch ein zweimaliges durchlaufen aller Elemente erfordern und wäre entsprechend ineffizient. Daher Definieren wir diesen äquivalent zum Additionsoperator. Beachten Sie, dass der zweite Parameter nie verändert werden muss und daher als Referenz auf const implementiert werden sollte.

template < typename T, typename Container >
inline matrix< T, Container >&
operator+=(matrix< T, Container >& lhs, matrix< T, Container > const& rhs){
    if(lhs.dimension() != rhs.dimension()){
        throw error::dimension_unequal("mitrax::operator+=<>()", lhs.dimension(), rhs.dimension());
    }
    typename matrix< T, Container >::iterator       target = lhs.begin();
    typename matrix< T, Container >::const_iterator add    = rhs.cbegin();
    for(;target != lhs.end(); ++target, ++add){
        *target += *add;
    }
    return lhs;
}

template < typename T, typename Container >
inline matrix< T, Container >&
operator-=(matrix< T, Container >& lhs, matrix< T, Container > const& rhs){
    if(lhs.dimension() != rhs.dimension()){
        throw error::dimension_unequal("mitrax::operator-=<>()", lhs.dimension(), rhs.dimension());
    }
    typename matrix< T, Container >::iterator       target = lhs.begin();
    typename matrix< T, Container >::const_iterator sub    = rhs.cbegin();
    for(;target != lhs.end(); ++target, ++sub){
        *target -= *sub;
    }
    return lhs;
}

template < typename T, typename Container >
inline matrix< T, Container > const
operator+(matrix< T, Container > lhs, matrix< T, Container > const& rhs){
    return lhs += rhs;
}

template < typename T, typename Container >
inline matrix< T, Container > const
operator-(matrix< T, Container > lhs, matrix< T, Container > const& rhs){
    return lhs -= rhs;
}

Skalarmultiplikation

Die Multiplikation einer Matrix mit einem Skalar, kann sowohl von rechts, als auch von links geschehen. Dabei werden alle Matrixelemente mit dem Skalar multipliziert. Die Division ist üblicherweise als Multiplikation mit dem Inversen definiert. Wenn also für zwei skalare Typen eine Division definiert ist, können wir auch eine Matrix elementweise durch einen Skalar dividieren. Die Division einer Matrix durch einen Skalar ergibt hingegen keinen sind.

Die Implementierung der Division über die Multiplikation mit dem Inversen des Skalars können wir, von den entstehenden Effizienzeinbußen einmal abgesehen, auch deshalb nicht vornehmen, weil wir für die Bildung des Inversen zumindest das neutrale Element für die Multiplikation vom Datentyp des Skalars kennen müssten. Falls Ihnen nicht bekannt ist, was das neutrale Element ist, es handelt sich dabei um den Wert, der bei einer Operation mit einem anderen Wert, diesen nicht verändert. Für ganze Zahlen ist dies zum Beispiel die Eins für die Multiplikation. Das neutrale Element der Addition ganzer Zahlen ist die Null.

Als Skalar könnten wir entweder nur Objekte vom Typ der Matrixelemente zulassen oder aber alle Typen die mit diesem Multipliziert werden können. Da letztes oft vorkommt, werden wir uns für diese Variante entscheiden. Denken Sie beispielsweise an eine Gleitkommamatrix, die mit einer Ganzzahl Multipliziert wird. Allerdings stehen wir an dieser Stelle vor einem Problem, denn wenn die Operatoren Templates sind, dann findet für die Parameter keine implizite Typumwandlung statt. Wenn wir aber einen Templateparameter für alle Skalartypen einführen, dann wird für jeden Skalartyp eine eigene Funktionsinstanz des Templates erzeugt und die implizite Umwandlung des Skalars wir für jedes Element neu ausgeführt. Eine gemeinsame Instantiierung mit dem Klassentemplate matrix können wir diesmal auch nicht realisieren, weil die Operatoren unabhängig von diesem eingebunden werden.

Was können wir also tun, damit die Umwandlung implizit bei der Übergabe des Skalarobjekts an den Operator erfolgt? Die Lösung besteht darin, die eigentliche Implementierung in zwei Funktionen im Namensraum detail auszulagern, welche nur den Matrixelementtyp als Skalartyp zulassen. Diese Funktionen werden dann von den verketteten Operatoren, unter Angabe der Templateparameter, aufgerufen. Auf diese Weise ist der Compiler nicht gezwungen, selbst nach der richtigen Funktionsinstanz der Templates zu suchen. Er benutzt einfach die Version, die wir von ihm verlangen. Da die Operatoren inline sind, dürfte der Compiler diese Aufrufe wegoptimieren. Falls mit der impliziten Typumwandlung ein Genauigkeitsverlust (z. B. float nach int) einhergeht, kann der Compiler immer noch warnen, da die Umwandlung weiterhin implizit ist.

namespace detail{
    template < typename T, typename Container >
    inline matrix< T, Container >&
    scalar_multiplication(matrix< T, Container >& lhs, T const& rhs){
        typedef typename matrix< T, Container >::iterator iterator;
        for(iterator target = lhs.begin(); target != lhs.end(); ++target){
            *target *= rhs;
        }
        return lhs;
    }

    template < typename T, typename Container >
    inline matrix< T, Container >&
    scalar_division(matrix< T, Container >& lhs, T const& rhs){
        typedef typename matrix< T, Container >::iterator iterator;
        for(iterator target = lhs.begin(); target != lhs.end(); ++target){
            *target /= rhs;
        }
        return lhs;
    }
}

template < typename T, typename Container, typename SkalarType >
inline matrix< T, Container >& operator*=(matrix< T, Container >& lhs, SkalarType const& rhs){
    return detail::scalar_multiplication< T, Container >(lhs, rhs);
}

template < typename T, typename Container, typename SkalarType >
inline matrix< T, Container >& operator/=(matrix< T, Container >& lhs, SkalarType const& rhs){
    return detail::scalar_division< T, Container >(lhs, rhs);
}

template < typename T, typename Container, typename SkalarType >
inline matrix< T, Container > const operator*(matrix< T, Container > lhs, SkalarType const& rhs){
    return lhs *= rhs;
}

template < typename T, typename Container, typename SkalarType >
matrix< T, Container > const operator*(SkalarType const& lhs, matrix< T, Container > rhs){
    return rhs *= lhs;
}

template < typename T, typename Container, typename SkalarType >
inline matrix< T, Container > const operator/(matrix< T, Container > lhs, SkalarType const& rhs){
    return lhs /= rhs;
}

Matrixmultiplikation

Bei der Matrixmultiplikation werden, wie der Name schon sagt, zwei Matrizen miteinander multipliziert. Zunächst nutzen wir typedef um die benötigten Typen unter einem kürzeren Namen verfügbar zu machen. Anschließend muss die Kompatibilität der Dimensionen der beiden Matrizen überprüft werden. Eine Matrixmultiplikation ist nur möglich, wenn die Anzahl der Spalten des linken Faktors und die Anzahl der Zeilen des rechten Faktors übereinstimmen. Ist dies nicht der Fall, werfen wir eine Ausnahme vom Typ error::dimension_incompatible.

Nun folgt die eigentliche Berechnung. Hierbei wird die Summe aus den Produkten der k'ten Elemente der i'ten Zeile des linken Faktors und der k'ten Elemente der j'ten Zeile des rechten Faktors Multipliziert, um das Element an der Position (i, j) zu erhalten. Wir legen uns also zunächst eine neue Matrix mit der richtigen Dimension an und berechnen anschließend die einzelnen Elemente. Das Ganze ist eine dreifach-Schleife. Das erste Element wird als Initialisierung genutzt und daher nicht in der innersten Schleife berechnet. Sofern Sie die verbale Beschreibung als zu schwierig empfanden, können Sie sich das Vorgehen auch anhand des Quelltextes deutlich machen.

template < typename T, typename Container >
inline matrix< T, Container >&
operator*=(matrix< T, Container >& lhs, matrix< T, Container > const& rhs){
    typedef typename matrix< T, Container >::size_type                    size_type;
    typedef typename matrix< T, Container >::row_const_proxy::iterator    row_iterator;
    typedef typename matrix< T, Container >::column_const_proxy::iterator column_iterator;

    // Dimensionskompatibilität prüfen
    if(lhs.columns() != rhs.rows()){
        throw error::dimension_incompatible(
            "mitrax::operator*() - Matrix multiplication requires lhs.columns() == rhs.rows()",
            lhs.dimension(), rhs.dimension()
        );
    }

    matrix< T, Container > result(lhs.rows(), rhs.columns());
    for(size_type column = size_type(); column < rhs.columns(); ++column){
        for(size_type row = size_type(); row < lhs.rows(); ++row){
            result[row][column] = lhs[row][0] * rhs[0][column];
            for(size_type line = size_type()+1; line < lhs.columns(); ++line){
                result[row][column] += lhs[row][line] * rhs[line][column];
            }
        }
    }

    return lhs = result;
}

template < typename T, typename Container >
inline matrix< T, Container > const
operator*(matrix< T, Container > lhs, matrix< T, Container > const& rhs){
    return lhs *= rhs;
}

Determinantenberechnung

Eine Determinante ist im mathematischen Sinne, eine Funktion, die einer quadratischen Matrix einen Skalar zuordnet. Fasst man die Matrix als Koeffizientenmatrix eines linearen Gleichungssystems auf, so ist dieses eindeutig lösbar, wenn die Determinante der Matrix Null ist. Für die Berechnung der Determinante einer Matrix gibt es sehr viele Lösungsverfahren. Einige sind sehr langsam, andere deutlich schneller. Wir werden uns natürlich tendenziell für ein schnelles Verfahren entscheiden. Dabei ist aber zu beachten, dass sich viele diese Verfahren die höhere Geschwindigkeit, mit einer geringeren numerischen Stabilität erkaufen. Die Rundungsfehler werden also deutlich größer. Daher ist es stark vom Anwendungszweck abhängig, ob man besser einen Algorithmus verwendet, der auf Geschwindigkeit optimiert oder einen, der numerisch besonders Stabil ist. Wenn der Nutzer mit einer der beiden Eigenschaften bei der Verwendung unserer Determinantenberechnung unzufrieden ist, dann kann und muss er sich eine eigene Version schreiben, die seinen jeweiligen Ansprüchen besser gerecht wird. Dabei kann er aber auf das komplette Repertoire von mitrax zurückgreifen und seine Erweiterung wird fast genauso zu verwenden sein, wie die nativen Funktionen. Wir werden dies in Kürze noch zeigen.

Unsere Implementierung wird nach dem gaußsches Eliminationsverfahren vorgehen. Dabei wird die Matrix in eine obere Dreiecksmatrix umgeformt. Für solche Matrizen ist die Determinante gleich dem Produkt der Diagonalelemente. Alle Element, die in einer oberen Dreiecksmatrix unterhalb beziehungsweise links der Diagonale stehen, sind Null. Die Umformung erfolgt durch mehrere Subtraktionen einer Zeile mit dem vielfachen einer anderen Zeile. Wird dies systematisch durchgeführt, so lassen sich für alle Elemente unterhalb der Diagonale Nullen erzeugen. Ist ein Diagonalelement selbst Null, so muss die entsprechende Zeile gegen eine Zeile getauscht werden, in der das Element in dieser Spalte nicht Null ist. Da wir nur ungern tatsächlich zwei komplette Zeilen austauschen möchten, werden wir den Zugriff mittels unserer Zeilenproxys realisieren und nötigenfalls einfach diese tauschen. Wenn keine Zeile gefunden werden kann, die an dieser Stelle ungleich Null ist, so ist die Determinante Null und die Berechnung kann beendet werden. Die Berechnungen für die Elemente, die im Verlauf der Berechnung zu Null gemacht werden, führen wir nicht aus. Wir wissen ja schon dass das Ergebnis Null sein wird und wir greifen auf die Elemente in der weitere Berechnung nicht noch einmal zu.

Für das Vertauschen von zwei Proxyobjekten werden wir einen Funktor verwenden, der sich intern gleich noch merkt, ob eine gerade oder ungerade Anzahl von Tauschvorgängen stattgefunden hat. Falls die Anzahl ungerade war, muss das Ergebnis am Ende nämlich noch einmal negiert werden. Zum Tauschen der Proxys setzen wir std::swap() ein. Wir verwenden dabei die using-Deklaration, auch wenn wir wissen, dass es keine spezielle Funktion dieses Namens für Proxyklassen gibt. Wenn wir dies konsequent immer tun, dann besteht nicht die Gefahr, dass wir es irgendwann einmal versehentlich anders machen. Da wir für diesen Algorithmus die Matrix selbst modifizieren, müssen wir natürlich auf einer Kopie der Übergeben arbeiten. Das kopieren erfolgt allerdings erst, nachdem wir sichergestellt haben, dass die Matrix auch quadratisch ist. Würden wir per „Call by value“ aufrufen und dann feststellen, dass ein Fehler geworfen werden muss, hätten wir die Kopie sonst nur erzeugt, um sie sogleich wieder zu Zerstören. Danach legen wir einen std::vector mit allen Zeilenproxys an. Der Zugriff ist dadurch, wie bei einem Matrixobjekt, durch zweimalige Anwendung des Indexoperators möglich, wir haben aber hier die Möglichkeit, die Proxys selbst, mittels std::swap, zu vertauschen. Die Implementierung sieht damit folgendermaßen aus:

namespace detail{
    template < typename Proxy >
    class proxy_swap_functor{
    public:
        proxy_swap_functor():odd_(false){}

        void operator()(Proxy& lhs, Proxy& rhs){
            using std::swap;
            swap(lhs, rhs);
            odd_ = !odd_;
        }

        operator bool(){return odd_;}

    private:
        bool odd_;
    };
}

template < typename T, typename Container >
inline T const det(matrix< T, Container > const& quad_matrix){
    typedef typename matrix< T, Container >::size_type                    size_type;
    typedef typename matrix< T, Container >::row_proxy                    row_proxy;
    typedef typename matrix< T, Container >::column_const_proxy::iterator column_iterator;

    if(quad_matrix.rows() != quad_matrix.columns()){
        throw error::not_quadratic("mitrax::det<>()", quad_matrix.dimension());
    }

    size_type size = quad_matrix.rows();
    matrix< T, Container > original(quad_matrix);

    // Liste von Zeilen-Proxys erstellen
    std::vector< row_proxy > lines;
    lines.reserve(size);
    for(size_type line = size_type(); line < size; ++line){
        lines.push_back(original.row(line));
    }

    detail::proxy_swap_functor< row_proxy > proxy_swap;

    // Obere Dreiecksmatrix bilden
    for(size_type line = size_type(); line < size; ++line){
        // Zeile mit einer folgenden tauschen, falls nötig
        if(lines[line][line] == identity< T >::additive()){
            bool status = proxy_swap;
            for(size_type row = line+1; row < size; ++row){
                if(lines[row][line] == identity< T >::additive()) continue;
                proxy_swap(lines[line], lines[row]);
                break;
            }
            // Kein tausch durchgeführt? ja -> det == 0
            if(status == proxy_swap){
                return identity< T >::additive();
            }
        }
        // Alle folgenden Zeilen in dieser Spalte 0 machen
        // Die Berechnung der Spalte selbst kann man sich sparen, deshalb line+1
        for(size_type row = line+1; row < size; ++row){
            T factor = (lines[row][line] / lines[line][line]);
            for(size_type column = line+1; column < size; ++column){
                lines[row][column] -= factor * lines[line][column];
            }
        }
    }

    // Diagonalelemente Multiplizieren
    T result(lines[size_type()][size_type()]);
    for(size_type line = size_type()+1; line < size; ++line){
        result *= lines[line][line];
    }

    // Negieren, falls Anzahl Vertauschungen ungerade
    if(proxy_swap) result = -result;

    return result;
}

Den Proxyvektor erstellen wir, indem wir zunächst die benötigten Menge an Speicher vorreservieren lassen und anschließen die Objekte mittels push_back() in diesem Speicher erstellen. Lassen Sie sich von dem Ausdruck identity< T >::additive() nicht verwirren, das wir sofort erläutert werden.

Neutrale Elemente

Wenn wir bisher über Nullen gesprochen haben, dann ist damit im Allgemeinen das neutrale Element bezüglich der Addition gemeint. Meist erzeugt der Standardkonstruktor ein solches Objekt, aber das muss nicht so sein. Das neutrale Element bezüglich der Multiplikation ist übrigens gleichzeitig auch jenes, das in einer Einheitsmatrix den Diagonalelementen entspricht. Bei diesen ist es noch viel unwahrscheinlicher, dass wir eine allgemeingültige Möglichkeit finden, um ein entsprechendes Objekt zu erzeugen. Sicher, für die eingebauten Datentypen und auch für unsere Bruchklasse oder für std::complex können wir ein Objekt standardkonstruiieren und mit dem int-Literal 1 addieren. Jedenfalls dann, wenn std::complex mit int Instantiiert wurde. In dem Moment, wo wir einen std::complex< double > haben, funktioniert dies schon nicht mehr. Noch schlimmer sieht es für unsere eigene Matrixklasse aus. Diese definiert eine entsprechende Addition nicht einmal in Sonderfällen, also auch nicht, wenn die Matrix nur ein Element vom Typ int enthält. Wir benötigen daher eine bessere Möglichkeit. Diese findet sich in Form sogenannter Traits.

Traits sind Helferklassen, die Informationen über Typen bereitstellen. In diesem Fall wollen wir ein Trait bereitstellen, das uns für seinen Templateparameter das neutrale Element bezüglich der Addition und bezüglich der Multiplikation bereitstellt. Der Nutzer kann dann eine Spezialisierung vornehmen, falls er einen Elementtyp hat, für den die Standardimplementierung versagt. Diese sieht im einfachsten Fall genau so aus, wie wir es eben angenommen haben. Das additive Neutrum entspricht einem standardkonstruiertem Objekt und das multiplikative Neutrum einem standardkonstruiertem Objekt plus den int-Literal 1. Diese kann implizit in alle anderen eingebauten Datentypen konvertiert werden und auch unsere Bruchklasse erlaubt mittels impliziter Typumwandlung eine entsprechende Addition. Beachten Sie, dass wir nicht einfach einen Konstruktor aufrufen sollen und diesem das Objekt 1 übergeben. Ein solcher Konsturktor muss sich nicht unbedingt auf einen Initialwert beziehen. Eine gute Klasse würde einen solchen Konstruktor zwar als explizit definieren, aber nicht alle Programmierer schreiben solch sauberen Code, um es mal vorsichtig auszudrücken. Die Traitklasse sieht wie folgt aus:

template < typename T >
struct identity{
    static T const& additive()      { static T const t;        return t; };
    static T const& multiplicative(){ static T const t(T()+1); return t; };
};

Das die Objekte innerhalb statischer Funktionen als statisch definiert sind, hat mehrere Gründe. Die Funktionen sind statisch, weil die Informationen nicht an ein bestimmtes Objekt gebunden sind. Sie geben Referenzen auf konstante Objekte zurück, weil nicht jedes mal das gleiche Objekt neu erzeugt werden soll. Die eigentlichen Objekte sind innerhalb der Funktionen und nicht innerhalb der Klasse als statisch deklariert, weil dadurch Erstens keine Unglücke bei der Initialisierungsreihenfolge passieren können und Zweitens die Initialisierung nur dann vorgenommen wird, wenn die entsprechende Funktion auch tatsächlich mindestens einmal im Programm aufgerufen wird. Auf die Geschichte mit der Initialisierungsreihenfolge wollen wir an dieser Stelle nicht näher eingehen, da dies numerische Datentypen ohnehin so gut wie nie betreffen dürfte. Wenn Sie allerdings keine Funktion haben, die die entsprechende Information für den Typ abfragt, dann wird auch niemals ein Objekt erzeugt. Der Speicher wird zum Programmstart alloziert, aber die Objekterzeugung erfolgt bei statischen Variablen, innerhalb von Funktionen, erst, beim Ersten Aufruf der Funktion. Speziell wenn Sie mit vielen verschiedenen Typen arbeiten, muss ihr Programm beim Start nicht erst viel Zeit damit verbringen, alles schon mal vorweg zu Berechnen, ohne dabei zu wissen, was später wirklich gebraucht wird.

Nun stehen wir aber immer noch vor dem Problem, dass die automatische Erzeugung einen multiplikationsneutralen Objekts für std::complex< double > nicht funktioniert. Wie jede gute numerische Templateklasse, stellt auch std::complex einen öffentlichen typedef namens value_type für den Datentyp seiner beiden Elemente bereit. Für diesen funktioniert die implizite Umwandlung. Den typedef selbst benötigen wir zwar nicht, aber er hat etwas mit der Lösung zu tun. Falls Sie bereits ein wenig Erfahrung im Umgang mit funktionaler Programmierung haben und vielleicht auch schon mit der Templatemetaprogrammierung in C++, dann sollte es Ihnen leicht fallen, die entsprechende Lösung anhand dieser Information zu finden. Wir werden das Problem natürlich gleich für alle Datentypen lösen, die als Templateparameter einen Datentyp für die internen Elemente erhalten.

Wir übergeben jetzt also neben einem Typ T, auch noch ein Template Numeric als Parameter, das seinerseits einen Typ erwartet. Dann spezialisieren wir unsere Traitklasse für Numeric< T >. Das additive Neutrum wird wieder durch den Standardkonstruktor erzeugt, das multiplikative erzeugen wir diesmal jedoch, durch eine Addition eines standardkonstruierten Objekts von Numeric< T > und dem multiplikativen Neutrum von T. Da ein entsprechender Typ eine solche implizite Umwandlung normalerweise anbietet oder einen entsprechenden Additionsoperator bereitstellt, entspricht das genau dem multiplikativen Neutrum von Numeric< T >.

template < typename T, template < typename U > class Numeric >
struct identity< Numeric< T > >{
    static Numeric< T > const& additive(){
        static Numeric< T > const t=T(); return t;
    };

    static Numeric< T > const& multiplicative(){
        static Numeric< T > const t(Numeric< T >() + identity< T >::multiplicative());
        return t;
    };
};

Damit haben wir auch beliebige Verbschachtelungen der genannten Typen abgedeckt. Nun bleibt aber noch die Frage zu klären, wie es mit den neutralen Elementen für unsere Matrix aussieht. Wenn wir schon einmal eine solche Traitklasse anbieten, wäre es natürlich schön, wenn diese auch mit unserem eigenen Datentyp funktionieren würde. Leider stoßen wir an dieser Stelle an eine Grenze, deren Überwindung schwierig ist. Die neutralen Elemente einer Matrix hängen von deren Dimension ab und diese wird erst zur Laufzeit festgelegt. Somit kann es also für Matrizen gleichen Typs, verschiedene neutrale Elemente geben, je nachdem, wie die Dimension des konkreten Objekt eben gerade aussieht. Wenn wir also eine automatische Erzeugung wollten, dann müssten wir immer ein Objekt, dass die passende Dimension bereitstellt mit an die Traitklasse übergeben. Das allein wäre noch nicht weiter schlimm, aber damit die Spezialisierung mit dem Original kompatibel bleibt, müssten auch für Datentypen Objekte übergeben werden, bei denen überhaupt keine Notwendigkeit dazu besteht. Beim Schreiben eines Template liegt schließlich noch keine Information darüber vor, ob der verwendete Typ compilierzeit- oder laufzeitabhängige neutrale Elemente besitzt.

Aus diesem Grund werden wir an dieser Stelle vorerst aufhören. Das Problem lässt sich zwar lösen, aber das würde schon deutlich weiter in die Templatemetaprogrammierung hineinreichen und ist daher nicht Thema dieses Abschnitts. Die mitrax-Operationen, die Informationen über neutrale Elemente abfragen, sind daher für verschachtelte Matrizen und ähnliche Typen nicht verfügbar.

Die beiden Klassen sollen in der Datei matrix.hpp stehen, da Sie auch benötigt werden könnten, wenn ein Benutzer sich entschließt, auf unsere Operatoren zu verzichten.

Berechnung der inversen Matrix

Die inverse einer Matrix ist jene Matrix, die mit der Ersten multipliziert, die Einheitsmatrix ergibt. Wenn Sie sich an die Matrixmultiplikation erinnern, wird bereits klar, dass eine invertierbare Matrix immer quadratisch sein muss. Außerdem muss die Determinante ungleich Null sein, damit eine Matrix erfolgreich invertiert werden kann. Die Berechnung ist dabei jener, bei der Determinante, sehr ähnlich, denn auch hier werden wir das gaußsche Eliminationsverfahren wieder einsetzen.

Genau wie bei der Determinantenberechnung bilden wir zunächst wieder eine obere Dreiecksmatrix. Diesmal werfen wir jedoch eine Ausnahme, wenn wir erkennen, dass die Determinante Null ist. Das ist jedoch noch nicht der einzige Unterschied. Bevor wir mit der eigentlichen Berechnung beginnen, legen wir eine Einheitsmatrix mit der Dimension unserer übergeben Matrix an und deren Zeilen formen wir vollkommen äquivalent zu denen, der übergeben Matrix, mit um. Die Vertauschung von Zeilen darf bei dieser Einheitsmatrix jedoch nicht durch die Vertauschung von Proxys durchgeführt werden. In diesem Fall müssen wir tatsächlich sämtliche Zeilenelemente gegeneinander Austauschen. Für diesen Zweck haben wir ja bereits die Funktion element_swap() geschrieben.

Nachdem wir eine obere Dreicksmatrix erreicht haben, müssen wir alle Zeilen durch das entsprechende Diagonalelement der übergebenen Matrix teilen, so dass alle Diagonalelemente dieser Matrix zu Eins werden. Das betrifft also auch die entsprechenden Zeilen der Einheitsmatrix. Anschließend müssen wir durch weitere Zeilensubtraktionen eine Diagonalmatrix bilden. Also eine Matrix, bei der nur die Diagonalelemente ungleich Null sind. Auch hierbei müssen wir die Subtraktionen wieder für die Zeilen beider Matrizen ausführen. Für die übergebene Matrix, müssen wir praktischerweise nicht wirklich Rechnen. Wie schon im ersten Schritt, wissen wir auch diesmal, dass die Elemente hinterher alle Null sein werden und dass wir später nicht noch einmal auf sie zugreifen. Nach Abschluss dieser Aktion, steht in der übergeben Matrix eine Einheitsmatrix und in der ehemaligen Einheitsmatrix die Inverse, des ursprünglich übergebenen matrix-Objekts. Da wir die unnötigen Berechnungen in der ersten Matrix nicht ausgeführt haben, steht in unserem Objekt keine echte Einheitsmatrix, aber wenn wir die Berechnungen ausgeführt hätten, dann Stände dort eine. Da unser Programm diese Objekt nun aber sowiso verwirft, und nur die Inverse zurückgibt, wäre es Unsinn, die Berechnungen wirklich bis zum entstehen einer Einheitsmatrix durchlaufen zu lassen.

template < typename T, typename Container >
inline matrix< T, Container > const inverse(matrix< T, Container > const& quad_matrix){
    typedef typename matrix< T, Container >::size_type                    size_type;
    typedef typename matrix< T, Container >::row_proxy                    row_proxy;
    typedef typename matrix< T, Container >::column_const_proxy::iterator column_iterator;

    if(quad_matrix.rows() != quad_matrix.columns()){
        throw error::not_quadratic("mitrax::inverse<>()", quad_matrix.dimension());
    }

    size_type size = quad_matrix.rows();
    matrix< T, Container > original(quad_matrix);

    // Einheitsmatrix erzeugen
    matrix< T, Container > result(quad_matrix.dimension());
    {
        typename matrix< T, Container >::iterator iter = result.begin();
        for(size_type i = size_type(); i < size; ++i){
            *iter = identity< T >::multiplicative();
            iter += size + 1;
        }
    }

    // Liste von Zeilen-Proxys erstellen
    std::vector< row_proxy > org_lines;
    org_lines.reserve(size);
    for(size_type line = size_type(); line < size; ++line){
        org_lines.push_back(original.row(line));
    }

    detail::proxy_swap_functor< row_proxy > proxy_swap;

    // Obere Dreiecksmatrix bilden
    for(size_type line = size_type(); line < size; ++line){
        // Zeile mit einer folgenden tauschen, falls nötig
        if(org_lines[line][line] == identity< T >::additive()){
            bool status = proxy_swap;
            for(size_type row = line+1; row < size; ++row){
                if(org_lines[row][line] == identity< T >::additive()) continue;
                proxy_swap(org_lines[line], org_lines[row]);
                element_swap(result[line], result[row]);
                break;
            }
            // Kein Tausch durchgeführt? ja -> det == 0 -> nicht invertierbar
            if(status == proxy_swap){
                throw error::singular_matrix("mitrax::inverse<>()");
            }
        }
        // Alle folgenden Zeilen in dieser Spalte 0 machen
        // (Elemente dieser Spalte müssen nicht berechnet werden)
        for(size_type row = line+1; row < size; ++row){
            T factor = (org_lines[row][line] / org_lines[line][line]);
            for(size_type column = line+1; column < size; ++column){
                org_lines[row][column] -= factor * org_lines[line][column];
            }
            for(size_type column = size_type(); column < size; ++column){
                result[row][column] -= factor * result[line][column];
            }
        }
        // Zeile normieren (Diagonalelemente müssen nicht berechnet werden)
        for(size_type column = line+1; column < size; ++column){
            org_lines[line][column] /= org_lines[line][line];
        }
        for(size_type column = size_type(); column < size; ++column){
            result[line][column] /= org_lines[line][line];
        }
    }

    // Einheitsmatrix bilden (Elemente in org_lines müssen nicht berechnet werden)
    for(size_type line = size-1; line > size_type(); --line){
        for(size_type row = size_type(); row < line; ++row){
            T factor(org_lines[row][line]);
            for(size_type column = size_type(); column < size; ++column){
                result[row][column] -= factor * result[line][column];
            }
        }
    }

    return result;
}

Wie sie sehen, eignet sich unser swap-Funktor in dieser Funktion immer noch bestens, um festzustellen, ob eine Zeilenvertauschung stattgefunden hat, auch wenn wir am Ende keine Auskunft über die Anzahl der Vertauschungen benötigen.

Vergleichsoperatoren

Für den Gleichheitsoperator werden wir zunächst prüfen ob die Dimensionen gleich sind. Wenn nicht, geben wir false zurück. Alternativ könnten wir an dieser Stelle auch eine Ausnahme werfen. Wenn die Dimensionen übereinstimmen, dann überprüfen wir elementweise auf Gleichheit. Stimmen alle Elemente überein, geben wir true zurück. Den Ungleichheitsoperator definieren wir auf die übliche Weise. Die übrigen vier Vergleichsoperatoren lassen sich für Matrizen nicht sinnvoll definieren, daher werden wir dies auch nicht tun.

template < typename T >
inline bool operator==(matrix< T, Container > const& lhs, matrix< T, Container > const& rhs){
    typedef typename matrix< T, Container >::const_iterator const_iterator;

    if(lhs.dimension() != rhs.dimension()){
        return false;
    }

    for(const_iterator lp = lhs.cbegin(), rp = rhs.cbegin(); lp != lhs.end(); ++lp, ++rp){
        if(*lp != *rp) return false;
    }

    return true;
}

template < typename T >
inline bool operator!=(dimension< T > const& lhs, dimension< T > const& rhs){
    return !(lhs == rhs);
}

Erweiterung um eigene Operationen

Wenn Sie mitrax um eigene Rechenoperation erweitern wollen, dann verwenden Sie einen eigenen Namensraum, in dem Ihre Erweiterungen stehen. Wir werden beispielhaft den Namensraum matrix_extension_username verwenden. Darin soll es eine Funktion transpose() geben, die ein mitrax::matrix-Objekt übernimmt, dieses in einer Kopie transponiert und das Ergebnis zurückgibt. Transponieren einer Matrix bedeutet, dass alle Zeilen zu Spalten und folglich alle Spalten zu Zeilen gemacht werden. Da es möglich ist, das diese Operation in einer späteren Version von mitrax nativ unterstützt wird und wir nicht wollen, dass es dadurch Compilierfehler gibt, darf die Erweiterung nicht im Namensraum mitrax stehen. Das gilt allgemein für alle Bibliotheken, Sie wissen schließlich nie, was noch kommt.

namespace matrix_extension_username{
    template < typename T, typename Container >
    inline mitrax::matrix< T, Container > const
    transpose(mitrax::matrix< T, Container > const& org){
        typedef typename mitrax::matrix< T, Container >::size_type size_type;

        mitrax::matrix< T, Container > result(org.columns(), org.rows());
        for(size_type row = size_type(); row < org.rows(); ++row){
            for(size_type column = size_type(); column < org.columns(); ++column){
                result[column][row] = org[row][column];
            }
        }

        return result;
    }
}

Diese Tatsache hat jedoch zufolge, dass die argumentabhängige Namenssuche nicht mehr funktioniert, weil die Funktion ja jetzt in einem anderen Namensraum steht, als ihr Parameter. Wenn Sie zwingend Ihre Version verwenden wollen, dann geben Sie Ihren Namensraum beim Aufruf von transpose() explizit an. Wenn Sie stattdessen lieber auf die native Version umsteigen wollen, sobald die Bibliothek diese unterstützt, dann verwenden Sie für Ihre Version eine using-Deklaration und rufen transpose() ohne Angabe des Bereichs auf. Wie bei der swap()-Funktion, wird der Compiler dann eine Funktion passenden Namens, die im Namensraum des Parameters deklariert ist, vorziehen.

// Benötigte includes

int main(){
    matrix m(3, 5);

    // Wertebelegung
    int j = 0;
    for(matrix::iterator i = m3.begin(); i != m3.end(); ++i){
        *i = j++;
    }

    // Aufruf der Erweiterung transpose() (falls nicht nativ vorhanden)
    using matrix_extension_username::transpose;
    std::cout << m3 << std::endl;
    std::cout << transpose(m3) << std::endl;
}

Falls Sie eine ganze Reihe von Erweiterungen dieser Art haben und alle nutzen möchten, sofern keine native Implementierung gleichen Namens existiert, dann können Sie auch eine using-Direktive verwenden, um alle Namen, aus Ihrem Namensraum für Erweiterungen verfügbar zu machen.