From ff8619c1c440fd1a7c8deaca763a33b8e1d53b5e Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 13 Mar 2022 21:32:47 +0200
Subject: [PATCH] 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.
---
 source/linal/dynamicmatrix.h | 30 +++++++++++++++---------------
 source/linal/matrix.h        | 24 ++++++++++++------------
 source/linal/matrixops.h     | 28 ++++++++++++++--------------
 3 files changed, 41 insertions(+), 41 deletions(-)

diff --git a/source/linal/dynamicmatrix.h b/source/linal/dynamicmatrix.h
index 1b191a1..de4511c 100644
--- 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 2711d40..bf21184 100644
--- 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 84c8012..0988acc 100644
--- 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;
 }
-- 
2.43.0