inverse does not work yet

src/matrix.hxx

@@ -177,26 +177,32 @@ template<typename TYPE, typename ARRAY=TYPE*> class MatrixBase {
     ///{
 
     TYPE det() {
-      // rule of sarrus
-      TYPE res(0);
-      assert(ROWS==COLUMNS); // not really necessary
-      // 1) terms to add
-      for (size_t i(0); i<COLUMNS; ++i) {
-        TYPE term(1);
-        for (size_t j(0); j<ROWS; ++j)
-          term *= at((i+j)%COLUMNS, j);
-        res += term;
-      }
-      // 2) terms to subtract
-      for (size_t i(0); i<COLUMNS; ++i) {
-        TYPE term(1);
-        for (size_t j(0); j<ROWS; ++j)
-          term *= at((i+j)%COLUMNS, ROWS-j-1);
-        res -= term;
-      }
+      TYPE res(gauss());
+      for (TYPE *p((TYPE*)(_c)+SIZE); --p>=(TYPE*)(_c); p-=COLUMNS) res *= *p;
       return res;
     }
 
+    TYPE gauss() {
+      /// calculate using gauss algorithmus
+      /// @see http://www.mathebibel.de/determinante-berechnen-nach-gauss
+      /// 1. normalize first line to first value
+      TYPE lambda(at(0, 0));
+      at(0, 0) = 1;
+      for (TYPE *p((TYPE*)(_c)+COLUMNS); p>(TYPE*)(_c)+1;) *--p/=lambda;
+      /// 2. nullify lower triangle
+      for (size_t column(0); column<COLUMNS-1; ++column) {
+        for (size_t row(column+1); row<ROWS; ++row) {
+          TYPE pivot(at(row, column));
+          if (pivot!=0) {
+            at(row, column) = 0;
+            for (size_t pos(column+1); pos<COLUMNS; ++pos)
+              at(row, pos) -= pivot*at(0, pos);
+          }
+        }
+      }
+      return lambda;
+    }
+
     ///}
     
     //................................................................methods
@@ -284,6 +290,42 @@ template<typename TYPE, size_t TROWS=0, size_t TCOLUMNS=0> class Matrix:
       return res;
     }
 
+    static Matrix i() {
+      Matrix res;
+      for (size_t row(0); row<TROWS&&row<TCOLUMNS; ++row) res(row, row) = 1;
+      return res;
+    }
+
+    Matrix& inv() {
+      /// calculate using gauss-jordan algorithmus
+      /// @see http://www.mathebibel.de/inverse-matrix-berechnen-nach-gauss-jordan
+      Matrix o(*this); // left side
+      *this = i(); // right side
+      /// 1. lower left part
+      for (size_t row(0); row<this->ROWS; ++row)
+        if (rows<this->COLUMNS) {
+          /// 2. normalize first line to first value
+          TYPE pivot(o(row, row));
+          if (pivot!=1) {
+            o(row, row) = 1;
+            for (size_t column(row+1); column<this->COLUMNS; ++column) o(row, column)/=pivot;
+            for (size_t column(row); column<this->COLUMNS; ++column) this->at(row, column)/=pivot;
+          }
+          /// 3. nullify lower triangle
+          for (size_t row(column+1); row<this->ROWS; ++row) {
+            TYPE pivot(o(row, column));
+            if (pivot!=0) {
+              o(row, column) = 0;
+              for (size_t pos(column+1); pos<this->COLUMNS; ++pos) o(row, pos) -= pivot*o(0, pos);
+              for (size_t pos(column); pos<this->COLUMNS; ++pos) this->at(row, pos) -= pivot*this->at(0, pos);
+            }
+          }
+        }
+      /// 4. nullify the upper triangle
+      
+      return *this;
+    }
+
     ///}
 };
 
@@ -306,6 +348,8 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
     
     Matrix(size_t rows, size_t columns):
         Parent(rows, columns) {
+      assert(rows>0);
+      assert(columns>0);
       Parent::_c = new TYPE[rows*columns];
       memset(Parent::_c, 0, Parent::MEM_SIZE);
     }
@@ -347,14 +391,6 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
       return Parent::operator=(o);
     }
 
-    Matrix t() const {
-      Matrix res(Parent::COLUMNS, Parent::ROWS);
-      for (size_t row(0); row<Parent::ROWS; ++row)
-        for (size_t column(0); column<Parent::COLUMNS; ++column)
-          res(column, row) = this->at(row, column);
-      return res;
-    }
-
     Matrix operator-() const {
       Matrix res(Parent::COLUMNS, Parent::ROWS);
       for (TYPE *to((TYPE*)(res._c)+Parent::SIZE); to>(TYPE*)(Parent::_c); --to)
@@ -374,16 +410,36 @@ template<typename TYPE> class Matrix<TYPE, 0, 0>: public MatrixBase<TYPE> {
 
     ///}
 
+    ///@name special operations
+    ///{
+    
+    Matrix t() const {
+      Matrix res(this->COLUMNS, this->ROWS);
+      for (size_t row(0); row<this->ROWS; ++row)
+        for (size_t column(0); column<this->COLUMNS; ++column)
+          res(column, row) = this->at(row, column);
+      return res;
+    }
+
+    Matrix i() const {
+      Matrix res(this->ROWS, this->COLUMNS);
+      for (size_t row(0); row<this->ROWS&&row<this->COLUMNS; ++row)
+        res(row, row) = 1;
+      return res;
+    }
+
+    ///}
+
     //................................................................methods
   protected:
 
     virtual void assert_check(const Matrix& o) const {
-      assert(o.ROWS==Parent::ROWS);
-      assert(o.COLUMNS==Parent::COLUMNS);
+      assert(o.ROWS==this->ROWS);
+      assert(o.COLUMNS==this->COLUMNS);
     }
 
     virtual bool check(const Matrix& o) const {
-      return o.ROWS==Parent::ROWS && o.COLUMNS==Parent::COLUMNS;
+      return o.ROWS==this->ROWS && o.COLUMNS==this->COLUMNS;
     }
 
     void copy_args(TYPE*) {}