Make gauss_jordan operate on columns instead of rows

[libs/math.git] / source / linal / dynamicmatrix.h

diff --git a/source/linal/dynamicmatrix.h b/source/linal/dynamicmatrix.h
index 1bbed275844c3fe06a60e4f09644c7c668bd4c5d..de4511c647e870eaf97376127c3195b8ccabe979 100644 (file)
--- a/source/linal/dynamicmatrix.h
+++ b/source/linal/dynamicmatrix.h
@@ -2,6 +2,7 @@
 #define MSP_LINAL_DYNAMICMATRIX_H_
 
 #include <algorithm>
+#include <ostream>
 #include "dynamicvector.h"
 #include "matrixops.h"
 
@@ -44,9 +45,9 @@ public:
        DynamicMatrix &operator+=(const DynamicMatrix &);
        DynamicMatrix &operator-=(const DynamicMatrix &);
 
-       DynamicMatrix &exchange_rows(unsigned, unsigned);
-       DynamicMatrix &multiply_row(unsigned, T);
-       DynamicMatrix &add_row(unsigned, unsigned, T);
+       DynamicMatrix &exchange_columns(unsigned, unsigned);
+       DynamicMatrix &multiply_column(unsigned, T);
+       DynamicMatrix &add_column(unsigned, unsigned, T);
 
        DynamicMatrix &invert();
 };
@@ -251,38 +252,38 @@ inline bool operator==(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
 }
 
 template<typename T>
-inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_rows(unsigned i, unsigned j)
+inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_columns(unsigned i, unsigned j)
 {
-       if(i>=rows_ || j>=rows_)
-               throw std::out_of_range("DynamicMatrix::exchange_rows");
+       if(i>=columns_ || j>=columns_)
+               throw std::out_of_range("DynamicMatrix::exchange_columns");
 
        using std::swap;
        for(unsigned k=0; k<columns_; ++k)
-               swap(element(i, k), element(j, k));
+               swap(element(k, i), element(k, j));
 
        return *this;
 }
 
 template<typename T>
-inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_row(unsigned i, T s)
+inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_column(unsigned i, T s)
 {
-       if(i>=rows_)
-               throw std::out_of_range("DynamicMatrix::multiply_row");
+       if(i>=columns_)
+               throw std::out_of_range("DynamicMatrix::multiply_column");
 
        for(unsigned k=0; k<columns_; ++k)
-               element(i, k) *= s;
+               element(k, i) *= s;
 
        return *this;
 }
 
 template<typename T>
-inline DynamicMatrix<T> &DynamicMatrix<T>::add_row(unsigned i, unsigned j, T s)
+inline DynamicMatrix<T> &DynamicMatrix<T>::add_column(unsigned i, unsigned j, T s)
 {
-       if(i>=rows_ || j>=rows_)
-               throw std::out_of_range("DynamicMatrix::exchange_rows");
+       if(i>=columns_ || j>=columns_)
+               throw std::out_of_range("DynamicMatrix::exchange_columns");
 
        for(unsigned k=0; k<columns_; ++k)
-               element(j, k) += element(i, k)*s;
+               element(k, j) += element(k, i)*s;
 
        return *this;
 }
@@ -307,7 +308,8 @@ inline DynamicMatrix<T> &DynamicMatrix<T>::invert()
        for(unsigned i=0; i<rows_; ++i)
                r(i, i) = T(1);
 
-       return gauss_jordan(*this, r);
+       gauss_jordan(*this, r);
+       return *this = r;
 }
 
 template<typename T>
@@ -317,6 +319,27 @@ inline DynamicMatrix<T> invert(const DynamicMatrix<T> &m)
        return r.invert();
 }
 
+template<typename T>
+inline std::ostream &operator<<(std::ostream &s, const DynamicMatrix<T> &m)
+{
+       s << "DynamicMatrix" << m.rows() << 'x' << m.columns() << '(';
+       for(unsigned i=0; i<m.columns(); ++i)
+       {
+               if(i)
+                       s << ", ";
+               s << '[';
+               for(unsigned j=0; j<m.rows(); ++j)
+               {
+                       if(j)
+                               s << ", ";
+                       s << m(j, i);
+               }
+               s << ']';
+       }
+       s << ')';
+       return s;
+}
+
 } // namespace LinAl
 } // namespace Msp