#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
