C++-Programmierung/ Eine Matrix-Bibliothek – mitrax/ simple operator.hpp

#ifndef _mitrax_mitrax_simple_operator_hpp_INCLUDED_
#define _mitrax_mitrax_simple_operator_hpp_INCLUDED_
/// \file simple_operator.hpp
///
/// \brief Operationen, die auf einer Matrix ausgeführt werden
///
/// Diese Datei stellt einige grundlegenden Operationen bereit, welche auf einer Matrix ausgeführt
/// werden können. Der Nutzer von mitrax kann in einem eigenen Namensraum beliebige weitere
/// Operationen schreiben. Wenn diese per <code>using</code>-Deklaration oder
/// <code>using</code>-Direktive verfügbar gemacht werden, lassen sie sich syntaktisch genau so
/// verwenden wie die hier definierten Funktionen. Ist eine Funktion sowohl als Erweiterung, als
/// auch als native (hier definierte Funktion) vorhanden, so wird der nativen Implementierung,
/// aufgrund der parameterabhängigen Namenssuchen, der Vorzug gegeben.

#include "exception.hpp"

namespace mitrax{

/// \brief Elementweise Negation einer Matrix
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;
}

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


/// \brief Additionszuweisung einer Matrix an eine andere
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;
}

/// \brief Subtraktionszuweisung einer Matrix an eine andere
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;
}


/// \brief Elementweise Addition zweier Matrizen
template < typename T, typename Container >
inline matrix< T, Container > const
operator+(matrix< T, Container > lhs, matrix< T, Container > const& rhs){
    return lhs += rhs;
}

/// \brief Elementweise Subtraktion zweier Matrizen
template < typename T, typename Container >
inline matrix< T, Container > const
operator-(matrix< T, Container > lhs, matrix< T, Container > const& rhs){
    return lhs -= rhs;
}


namespace detail{

/// \brief Multiplikation aller Elemente einer Matrix mit einem Skalar (Helferfunktion)
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;
}

/// \brief Division aller Elemente einer Matrix durch einen Skalar (Helferfunktion)
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;
}

}

/// \brief Multiplikation aller Elemente einer Matrix mit einem Skalar
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);
}

/// \brief Division aller Elemente einer Matrix durch einen Skalar
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);
}

/// \brief Matrixmultiplikation mit Zuweisung
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;
}

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

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

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

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

namespace detail{

/// \brief Funktor zum Vertauschen zweier Matrixproxys
template < typename Proxy >
class proxy_swap_functor{
public:
    /// \brief Initialisierung mit <code>false</code>
    proxy_swap_functor():odd_(false){}

    /// \brief Vertauscht zwei Proxys und wechselt intern seinen Wert
    void operator()(Proxy& lhs, Proxy& rhs){
        using std::swap;
        swap(lhs, rhs);
        odd_ = !odd_;
    }

    /// \brief <code>true</code> für eine ungerade Anzahl von erfolgten Vertauschungen
    operator bool(){return odd_;}

private:
    /// \brief <code>true</code> für eine ungerade Anzahl von erfolgten Vertauschungen
    bool odd_;
};

}


/// \brief Determinantenberechnung (nur quadratische Matrizen)
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;
}

/// \brief Inverse Matrix (nur reguläre quadratische Matrizen)
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;
}


/// \brief <code>true</code> für gleiche Dimension und Elemente
template < typename T, typename Container >
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;
}

/// \brief <code>true</code> für ungleiche Dimension oder ungleiche Elemente
template < typename T, typename Container >
inline bool operator!=(dimension< T > const& lhs, dimension< T > const& rhs){
    return !(lhs == rhs);
}


}

#endif