less use of transpose

include/Gauss.h

@@ -10,11 +10,19 @@ class Matrix;
 namespace Gauss {
 
 /**
- * \brief Echelonne une matrice en utilisant l'algorithme de Gauss-Jordan
- * \param mat La matrice à échelonner
- * \param reduite Mets des 0 au dessus des pivots
- * \param normalise Mets les pivots à 1
+ * \brief Echelonne une matrice en ligne en utilisant l'algorithme de Gauss-Jordan
+ * \param a_Matrix La matrice à échelonner
+ * \param a_Reduite Mets des 0 au dessus des pivots
+ * \param a_Normalise Mets les pivots à 1
  */
 void GaussJordan(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise);
 
+/**
+ * \brief Echelonne une matrice en colonne en utilisant l'algorithme de Gauss-Jordan
+ * \param a_Matrix La matrice à échelonner
+ * \param a_Reduite Mets des 0 au dessus des pivots
+ * \param a_Normalise Mets les pivots à 1
+ */
+void GaussJordanColumn(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise);
+
 }  // namespace Gauss

include/Matrix.h

@@ -66,12 +66,19 @@ class Matrix {
 	void Transpose();
 
 	/**
-	 * \brief Augmente la matrice actuelle avec une autre
+	 * \brief Augmente la matrice actuelle à droite avec une autre
 	 * \param a_Right Une matrice avec le bon nombre de lignes
 	 * \pre Les deux matrices doivent avoir le même nombre de lignes
 	 */
 	void Augment(const Matrix& a_Right);
 
+	/**
+	 * \brief Augmente la matrice actuelle en dessous avec une autre
+	 * \param a_Bottom Une matrice avec le bon nombre de colonnes
+	 * \pre Les deux matrices doivent avoir le même nombre de colonnes
+	 */
+	void AugmentBottom(const Matrix& a_Bottom);
+
 	/**
 	 * \brief Affecte tous les coefficients de la matrice à un élément
 	 * \param a_Element L'élément à affecter

include/Solver.h

@@ -18,14 +18,14 @@ class Solver {
 	 * \param a_Matrix La matrice à traiter
 	 * \return L'espace vectoriel correspondant
 	 */
-	Vect Image(const Matrix& a_Matrix) const;
+	Vect Image(Matrix&& a_Matrix) const;
 
 	/**
 	 * \brief Calcule le noyau d'une matrice
 	 * \param a_Matrix La matrice à traiter
 	 * \return L'espace vectoriel correspondant
 	 */
-	Vect Kernel(const Matrix& a_Matrix) const;
+	Vect Kernel(Matrix&& a_Matrix) const;
 
 	/**
 	 * \brief Résout le système rectangulaire de la forme AX=B, avec X et B, des vecteurs colonne.
@@ -33,12 +33,12 @@ class Solver {
 	 * \param a_VectorB La matrice colonne jouant le rôle de B
 	 * \return L'espace affine associé
 	 */
-	VectAffine RectangularSystem(const Matrix& a_MatrixA, const Matrix& a_VectorB) const;
+	VectAffine RectangularSystem(Matrix&& a_MatrixA, const Matrix& a_VectorB) const;
 
 	/**
 	 * \brief Calcule le rang d'une matrice
 	 * \param a_Matrix La matrice à traiter
 	 * \note Ceci équivaut à \code Image(a_Matrix).GetCardinal() \endcode
 	 */
-	std::size_t Rank(const Matrix& a_Matrix) const;
+	std::size_t Rank(Matrix&& a_Matrix) const;
 };

src/Gauss.cpp

@@ -10,12 +10,24 @@ static void SwapLines(Matrix& mat, std::size_t line1, std::size_t line2) {
 	}
 }
 
+static void SwapColumns(Matrix& mat, std::size_t column1, std::size_t column2) {
+	for (std::size_t k = 0; k < mat.GetRawCount(); k++) {
+		std::swap(mat.at(k, column1), mat.at(k, column2));
+	}
+}
+
 static void DivideLine(Matrix& mat, std::size_t line, Matrix::Element number) {
 	for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {
 		mat.at(line, j) /= number;
 	}
 }
 
+static void DivideColumn(Matrix& mat, std::size_t column, Matrix::Element number) {
+	for (std::size_t j = 0; j < mat.GetRawCount(); j++) {
+		mat.at(j, column) /= number;
+	}
+}
+
 static int FirstNotNullElementIndexOnColumn(Matrix& mat, std::size_t column, std::size_t startLine = 0) {
 	for (std::size_t i = startLine; i < mat.GetRawCount(); i++) {
 		if (!IsEqualZero(mat.at(i, column))) {
@@ -25,6 +37,15 @@ static int FirstNotNullElementIndexOnColumn(Matrix& mat, std::size_t column, std
 	return -1;
 }
 
+static int FirstNotNullElementIndexOnLine(Matrix& mat, std::size_t line, std::size_t startColumn = 0) {
+	for (std::size_t i = startColumn; i < mat.GetColumnCount(); i++) {
+		if (!IsEqualZero(mat.at(line, i))) {
+			return i;
+		}
+	}
+	return -1;
+}
+
 static void SimplifyLine(Matrix& mat, std::size_t line, std::size_t pivot_line, std::size_t pivot_column) {
 	const Matrix::Element pivot = mat.at(pivot_line, pivot_column);
 	const Matrix::Element anul = mat.at(line, pivot_column);
@@ -34,6 +55,15 @@ static void SimplifyLine(Matrix& mat, std::size_t line, std::size_t pivot_line,
 	}
 }
 
+static void SimplifyColumn(Matrix& mat, std::size_t column, std::size_t pivot_line, std::size_t pivot_column) {
+	const Matrix::Element pivot = mat.at(pivot_line, pivot_column);
+	const Matrix::Element anul = mat.at(pivot_line, column);
+
+	for (std::size_t j = 0; j < mat.GetRawCount(); j++) {
+		mat.at(j, column) = mat.at(j, column) * pivot - mat.at(j, pivot_column) * anul;
+	}
+}
+
 void GaussJordan(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise) {
 	int indice_ligne_pivot = -1;
 
@@ -66,4 +96,36 @@ void GaussJordan(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise) {
 	}
 }
 
+void GaussJordanColumn(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise) {
+	int indice_colonne_pivot = -1;
+
+	for (std::size_t j = 0; j < a_Matrix.GetRawCount(); j++) {
+
+		int indice_colonne_pivot_trouve = FirstNotNullElementIndexOnLine(a_Matrix, j, indice_colonne_pivot + 1);
+
+		if (indice_colonne_pivot_trouve < 0)  // ligne de 0
+			continue;						  // on regarde la prochaine ligne
+
+		indice_colonne_pivot++;
+
+		if (indice_colonne_pivot_trouve != indice_colonne_pivot) {
+			SwapColumns(a_Matrix, indice_colonne_pivot_trouve, indice_colonne_pivot);
+		}
+
+		Matrix::Element pivot = a_Matrix.at(j, indice_colonne_pivot);
+
+		if (a_Normalise) {
+			DivideColumn(a_Matrix, indice_colonne_pivot, pivot);
+		}
+
+		// On simplifie les autres lignes
+		for (std::size_t i = (a_Reduite ? 0 : j); i < a_Matrix.GetColumnCount(); i++) {
+			// Pour les lignes autre que la ligne pivot
+			if (i != static_cast<std::size_t>(indice_colonne_pivot)) {
+				SimplifyColumn(a_Matrix, i, j, indice_colonne_pivot);
+			}
+		}
+	}
+}
+
 }  // namespace Gauss

src/Matrix.cpp

@@ -97,6 +97,25 @@ void Matrix::Augment(const Matrix& a_Right) {
 	*this = temp;
 }
 
+void Matrix::AugmentBottom(const Matrix& a_Bottom) {
+	assert(a_Bottom.m_Columns == m_Columns);
+	Matrix temp {m_Raws + a_Bottom.GetRawCount(), m_Columns};
+
+	for (std::size_t i = 0; i < m_Raws; i++) {
+		for (std::size_t j = 0; j < m_Columns; j++) {
+			temp.at(i, j) = at(i, j);
+		}
+	}
+
+	for (std::size_t i = 0; i < a_Bottom.GetRawCount(); i++) {
+		for (std::size_t j = 0; j < GetColumnCount(); j++) {
+			temp.at(i + GetRawCount(), j) = a_Bottom.at(i, j);
+		}
+	}
+
+	*this = temp;
+}
+
 Matrix Matrix::operator+(const Matrix& a_Other) const {
 	assert(GetColumnCount() == a_Other.GetColumnCount() && GetRawCount() == a_Other.GetRawCount());
 

src/Solver.cpp

@@ -2,37 +2,34 @@
 
 #include "Gauss.h"
 
-Vect Solver::Image(const Matrix& a_Matrix) const {
-	Matrix result = a_Matrix;
-	result.Transpose();
-	Gauss::GaussJordan(result, true, true);
-	result.Transpose();
-	return {result};
+Vect Solver::Image(Matrix&& a_Matrix) const {
+	Gauss::GaussJordanColumn(a_Matrix, true, true);
+	return {a_Matrix};
 }
 
 // https://en.wikipedia.org/wiki/Kernel_(linear_algebra)#Computation_by_Gaussian_elimination
-Vect Solver::Kernel(const Matrix& a_Matrix) const {
-	Matrix result = a_Matrix;
-	result.Transpose();
-	result.Augment(Matrix::Identity(result.GetRawCount()));
-	Gauss::GaussJordan(result, true, true);
-	result.Transpose();
+Vect Solver::Kernel(Matrix&& a_Matrix) const {
+	std::size_t matrixRawCount = a_Matrix.GetRawCount();
+	std::size_t matrixColumnCount = a_Matrix.GetColumnCount();
+
+	a_Matrix.AugmentBottom(Matrix::Identity(a_Matrix.GetColumnCount()));
+	Gauss::GaussJordanColumn(a_Matrix, true, true);
 
 	// nombre de colonnes non nulles
-	std::size_t origine_colonne = Vect(result.SubMatrix(0, 0, a_Matrix.GetRawCount(), a_Matrix.GetColumnCount())).GetCardinal();
+	std::size_t origine_colonne = Vect(a_Matrix.SubMatrix(0, 0, matrixRawCount, matrixColumnCount)).GetCardinal();
 
-	return {result.SubMatrix(a_Matrix.GetRawCount(), origine_colonne, result.GetRawCount() - a_Matrix.GetRawCount(),
-		result.GetColumnCount() - origine_colonne)};
+	return {a_Matrix.SubMatrix(
+		matrixRawCount, origine_colonne, a_Matrix.GetRawCount() - matrixRawCount, a_Matrix.GetColumnCount() - origine_colonne)};
 }
 
-VectAffine Solver::RectangularSystem(const Matrix& a_MatrixA, const Matrix& a_VectorB) const {
+VectAffine Solver::RectangularSystem(Matrix&& a_MatrixA, const Matrix& a_VectorB) const {
 	Matrix mat = a_MatrixA;
 	mat.Augment(a_VectorB);
 	Gauss::GaussJordan(mat, true, true);
 
 	Solver solver;
 
-	Vect noyau = solver.Kernel(a_MatrixA);
+	Vect noyau = solver.Kernel(std::move(a_MatrixA));
 	Matrix origin = mat.SubMatrix(0, mat.GetColumnCount() - 1, mat.GetRawCount(), 1);
 
 	// on rajoute des 0 si il faut
@@ -49,6 +46,6 @@ VectAffine Solver::RectangularSystem(const Matrix& a_MatrixA, const Matrix& a_Ve
 	return {noyau, fullOrigin};
 }
 
-std::size_t Solver::Rank(const Matrix& a_Matrix) const {
-	return Image(a_Matrix).GetCardinal();
+std::size_t Solver::Rank(Matrix&& a_Matrix) const {
+	return Image(std::move(a_Matrix)).GetCardinal();
 }

src/Vect.cpp

@@ -75,9 +75,9 @@ Matrix Vect::GetLinearSystem() const {
 	vect.Transpose();
 
 	Solver solver;
-	vect = solver.Kernel(vect).m_Data;
-	vect.Transpose();
-	return vect;
+	Matrix result = solver.Kernel(std::move(vect)).m_Data;
+	result.Transpose();
+	return result;
 }
 
 std::size_t Vect::GetDimension() const {

src/main.cpp

@@ -25,8 +25,8 @@ void test() {
 
 	Solver solver;
 
-	Vect image = solver.Image(mat2);
-	Vect noyau = solver.Kernel(mat2);
+	Vect image = solver.Image(Matrix{mat2});
+	Vect noyau = solver.Kernel(Matrix{mat2});
 
 	std::cout << "\tImage :\n";
 	Print(image);

test/test_assert.h

@@ -2,21 +2,20 @@
 
 /**
  * \file Test.h
- * \brief File containing unit testing utilities
+ * \brief Contient une assertion utilisable avec les optimisations
  */
 
-#include <cstdlib>
 #include <iostream>
 
 /**
  * \def TEST_SUCCESSFUL
- * \brief Used in tests to indicate that a test was successful
+ * \brief Indique que le test a été passé
  */
 #define TEST_SUCCESSFUL 0
 
 /**
  * \def TEST_FAILED
- * \brief Used in tests to indicate that a test failed
+ * \brief Indique que le test a échoué
  */
 #define TEST_FAILED 1
 
@@ -30,9 +29,9 @@
 
 /**
  * \def test_assert
- * \param ... The expression to evaluate
- * \brief Evaluates the expression and exits the program if not valid.
- * \note It works like a basic assert() but also in release mode
+ * \param ... L'expression à évaluer
+ * \brief Evalue une expression et arrête le programme si elle n'est pas valide
+ * \note Cette macro équivaut à assert() mais fonctionne également avec les optimisations activées
  */
 #define test_assert(...)                                                                                                              \
 	if (!static_cast<bool>(__VA_ARGS__)) {                                                                                            \

test/test_random_kernel.cpp

@@ -37,7 +37,9 @@ static void Test() {
 		}
 	}
 
-	Vect kernel = solver.Kernel(matrix);
+	Matrix copy = matrix;
+
+	Vect kernel = solver.Kernel(std::move(copy));
 
 	Matrix nullVector {matrix.GetRawCount(), 1};
 	nullVector.Fill(0.0);

test/test_solver.cpp

@@ -18,7 +18,7 @@ void TestRectangular() {
 
 	Solver solver;
 
-	std::cout << solver.RectangularSystem(mat2, Matrix::ColumnVector({1, 2})).GetLinearSystem() << std::endl;
+	std::cout << solver.RectangularSystem(std::move(mat2), Matrix::ColumnVector({1, 2})).GetLinearSystem() << std::endl;
 	std::cout << aff.GetLinearSystem() << std::endl;
 }
 
@@ -39,8 +39,10 @@ void TestKernelImage() {
 
 		Solver solver;
 
-		test_assert(solver.Image(mat) == image);
-		test_assert(solver.Kernel(mat) == noyau);
+		Matrix copy = mat;
+
+		test_assert(solver.Image(std::move(copy)) == image);
+		test_assert(solver.Kernel(std::move(mat)) == noyau);
 	}
 }