Make gauss_jordan operate on columns instead of rows

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

diff --git a/source/linal/matrixops.h b/source/linal/matrixops.h
index 84c8012dd49addb539ee5d5f5ebcc4dfaf745fe1..0988accd6d44798ab906851fc180acf4455cba3e 100644 (file)
--- a/source/linal/matrixops.h
+++ b/source/linal/matrixops.h
@@ -22,37 +22,37 @@ inline T &gauss_jordan(T &m, T &r)
        typedef typename T::ElementType V;
        using std::abs;
 
-       for(unsigned i=0; i<m.rows(); ++i)
+       for(unsigned i=0; i<m.columns(); ++i)
        {
                unsigned pivot = i;
-               for(unsigned j=i+1; j<m.rows(); ++j)
-                       if(abs(m.element(j, i))>abs(m.element(pivot, i)))
+               for(unsigned j=i+1; j<m.columns(); ++j)
+                       if(abs(m.element(i, j))>abs(m.element(i, pivot)))
                                pivot = j;
 
-               if(m.element(pivot, i)==V(0))
+               if(m.element(i, pivot)==V(0))
                        throw not_invertible();
 
                if(pivot!=i)
                {
-                       m.exchange_rows(i, pivot);
-                       r.exchange_rows(i, pivot);
+                       m.exchange_columns(i, pivot);
+                       r.exchange_columns(i, pivot);
                }
 
-               for(unsigned j=i+1; j<m.rows(); ++j)
+               for(unsigned j=i+1; j<m.columns(); ++j)
                {
-                       V a = -m.element(j, i)/m.element(i, i);
-                       m.add_row(i, j, a);
-                       r.add_row(i, j, a);
+                       V a = -m.element(i, j)/m.element(i, i);
+                       m.add_column(i, j, a);
+                       r.add_column(i, j, a);
                }
 
                V a = V(1)/m.element(i, i);
-               m.multiply_row(i, a);
-               r.multiply_row(i, a);
+               m.multiply_column(i, a);
+               r.multiply_column(i, a);
        }
 
-       for(unsigned i=m.rows(); i-->0; )
+       for(unsigned i=m.columns(); i-->0; )
                for(unsigned j=i; j-->0; )
-                       r.add_row(i, j, -m.element(j, i));
+                       r.add_column(i, j, -m.element(i, j));
 
        return r;
 }