Make gauss_jordan operate on columns instead of rows

Since matrices are stored in column-major order, this makes the algorithm
more amenable to possible simdification.

diff --git a/source/linal/dynamicmatrix.h b/source/linal/dynamicmatrix.h
index 1b191a1b3f9a5d4696c4d535e94553e6b96f93e6..de4511c647e870eaf97376127c3195b8ccabe979 100644 (file)
--- a/source/linal/dynamicmatrix.h
+++ b/source/linal/dynamicmatrix.h
@@ -45,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();
 };
@@ -252,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;
 }

diff --git a/source/linal/matrix.h b/source/linal/matrix.h
index 2711d40a637489c1bfa3cc9676c2fb8043f978d4..bf21184092046c71bc1973ae1575deff51406fb7 100644 (file)
--- a/source/linal/matrix.h
+++ b/source/linal/matrix.h
@@ -51,9 +51,9 @@ public:
        Matrix &operator+=(const Matrix &);
        Matrix &operator-=(const Matrix &);
 
-       Matrix &exchange_rows(unsigned, unsigned);
-       Matrix &multiply_row(unsigned, T);
-       Matrix &add_row(unsigned, unsigned, T);
+       Matrix &exchange_columns(unsigned, unsigned);
+       Matrix &multiply_column(unsigned, T);
+       Matrix &add_column(unsigned, unsigned, T);
 };
 
 template<typename T, unsigned M, unsigned N>
@@ -227,27 +227,27 @@ inline bool operator==(const Matrix<T, M, N> &a, const Matrix<T, M, N> &b)
 }
 
 template<typename T, unsigned M, unsigned N>
-inline Matrix<T, M, N> &Matrix<T, M, N>::exchange_rows(unsigned i, unsigned j)
+inline Matrix<T, M, N> &Matrix<T, M, N>::exchange_columns(unsigned i, unsigned j)
 {
        using std::swap;
-       for(unsigned k=0; k<N; ++k)
-               swap(element(i, k), element(j, k));
+       for(unsigned k=0; k<M; ++k)
+               swap(element(k, i), element(k, j));
        return *this;
 }
 
 template<typename T, unsigned M, unsigned N>
-inline Matrix<T, M, N> &Matrix<T, M, N>::multiply_row(unsigned i, T s)
+inline Matrix<T, M, N> &Matrix<T, M, N>::multiply_column(unsigned i, T s)
 {
-       for(unsigned k=0; k<N; ++k)
-               element(i, k) *= s;
+       for(unsigned k=0; k<M; ++k)
+               element(k, i) *= s;
        return *this;
 }
 
 template<typename T, unsigned M, unsigned N>
-inline Matrix<T, M, N> &Matrix<T, M, N>::add_row(unsigned i, unsigned j, T s)
+inline Matrix<T, M, N> &Matrix<T, M, N>::add_column(unsigned i, unsigned j, T s)
 {
-       for(unsigned k=0; k<N; ++k)
-               element(j, k) += element(i, k)*s;
+       for(unsigned k=0; k<M; ++k)
+               element(k, j) += element(k, i)*s;
        return *this;
 }
 

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;
 }