more optimizations

include/Gauss.h

@@ -17,12 +17,4 @@ namespace Gauss {
  */
 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

@@ -17,6 +17,7 @@
 class Matrix {
   public:
 	typedef long double Element;
+	typedef std::vector<Element>::iterator iterator;
 
   private:
 	std::size_t m_Raws;
@@ -145,6 +146,11 @@ class Matrix {
 	 * construit une matrice de 1 ligne et 4 colonnes de coordonnées (1, 2, 3, 4)
 	 */
 	static Matrix RawVector(std::initializer_list<Element>&& a_Elements);
+
+	iterator begin();
+	iterator end();
+	
+	iterator GetLineIterator(std::size_t a_Raw);
 };
 
 template <typename T>

include/Vect.h

@@ -21,7 +21,7 @@ class Vect {
 	 * Les colonnes de 0 sont ignorées
 	 * \param a_Matrix Une matrice échelonnée.
 	 */
-	Vect(const Matrix& a_Matrix);
+	Vect(Matrix&& a_Matrix);
 
 	/**
 	 * \brief Permet d'obtenir le ieme vecteur de la base

src/Gauss.cpp

@@ -2,30 +2,20 @@
 
 #include "Matrix.h"
 
+#include <algorithm>
+#include <execution>
+#include <ranges>
+
 namespace Gauss {
 
 static void SwapLines(Matrix& mat, std::size_t line1, std::size_t line2) {
-	for (std::size_t k = 0; k < mat.GetColumnCount(); k++) {
-		std::swap(mat.at(line1, k), mat.at(line2, k));
-	}
-}
-
-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));
-	}
+	std::swap_ranges(
+		std::execution::par_unseq, mat.GetLineIterator(line1), mat.GetLineIterator(line1 + 1), mat.GetLineIterator(line2));
 }
 
 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;
-	}
+	std::transform(std::execution::par_unseq, mat.GetLineIterator(line), mat.GetLineIterator(line + 1), mat.GetLineIterator(line),
+		[number](Matrix::Element e) { return e /= number; });
 }
 
 static int FirstNotNullElementIndexOnColumn(Matrix& mat, std::size_t column, std::size_t startLine = 0) {
@@ -37,34 +27,18 @@ 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);
 
-	for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {
+	auto range = std::views::iota(static_cast<std::size_t>(0), mat.GetColumnCount());
+
+	std::for_each(std::execution::par_unseq, range.begin(), range.end(), [&mat, pivot, anul, line, pivot_line](std::size_t j) {
 		mat.at(line, j) = mat.at(line, j) * pivot - mat.at(pivot_line, j) * anul;
-	}
+	});
 }
 
-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) {
+static void GaussJordanReduced(Matrix& a_Matrix, bool a_Normalise) {
 	int indice_ligne_pivot = -1;
 
 	for (std::size_t j = 0; j < a_Matrix.GetColumnCount(); j++) {
@@ -86,46 +60,54 @@ void GaussJordan(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise) {
 			DivideLine(a_Matrix, indice_ligne_pivot, pivot);
 		}
 
+		auto range = std::views::iota(static_cast<std::size_t>(0), a_Matrix.GetRawCount());
+
 		// On simplifie les autres lignes
-		for (std::size_t i = (a_Reduite ? 0 : j); i < a_Matrix.GetRawCount(); i++) {
-			// Pour les lignes autre que la ligne pivot
+		std::for_each(std::execution::par_unseq, range.begin(), range.end(), [&a_Matrix, j, indice_ligne_pivot](std::size_t i) {
 			if (i != static_cast<std::size_t>(indice_ligne_pivot)) {
 				SimplifyLine(a_Matrix, i, indice_ligne_pivot, j);
 			}
+		});
 	}
 }
+
+static void GaussJordanTriangular(Matrix& a_Matrix, bool a_Normalise) {
+	int indice_ligne_pivot = -1;
+
+	for (std::size_t j = 0; j < a_Matrix.GetColumnCount(); j++) {
+
+		int indice_ligne_pivot_trouve = FirstNotNullElementIndexOnColumn(a_Matrix, j, indice_ligne_pivot + 1);
+
+		if (indice_ligne_pivot_trouve < 0)	// colonne de 0
+			continue;						// on regarde la prochaine colonne
+
+		indice_ligne_pivot++;
+
+		if (indice_ligne_pivot_trouve != indice_ligne_pivot) {
+			SwapLines(a_Matrix, indice_ligne_pivot_trouve, indice_ligne_pivot);
 		}
 
-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);
+		Matrix::Element pivot = a_Matrix.at(indice_ligne_pivot, j);
 
 		if (a_Normalise) {
-			DivideColumn(a_Matrix, indice_colonne_pivot, pivot);
+			DivideLine(a_Matrix, indice_ligne_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);
-			}
+		auto range = std::views::iota(static_cast<std::size_t>(indice_ligne_pivot + 1), a_Matrix.GetRawCount());
+
+		// On simplifie les autres lignes après la ligne du pivot
+		std::for_each(std::execution::par_unseq, range.begin(), range.end(),
+			[&a_Matrix, indice_ligne_pivot, j](std::size_t i) { 
+				SimplifyLine(a_Matrix, i, indice_ligne_pivot, j); 
+			});
 	}
 }
+
+void GaussJordan(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise) {
+	if (a_Reduite)
+		GaussJordanReduced(a_Matrix, a_Normalise);
+	else
+		GaussJordanTriangular(a_Matrix, a_Normalise);
 }
 
 }  // namespace Gauss

src/Matrix.cpp

@@ -187,3 +187,15 @@ Matrix Matrix::SubMatrix(
 
 	return result;
 }
+
+Matrix::iterator Matrix::begin() {
+	return m_Data.begin();
+}
+
+Matrix::iterator Matrix::end() {
+	return m_Data.end();
+}
+
+Matrix::iterator Matrix::GetLineIterator(std::size_t a_Raw) {
+	return m_Data.begin() + a_Raw * GetColumnCount();
+}

src/Solver.cpp

@@ -3,8 +3,10 @@
 #include "Gauss.h"
 
 Vect Solver::Image(Matrix&& a_Matrix) const {
-	Gauss::GaussJordanColumn(a_Matrix, true, true);
-	return {a_Matrix};
+	a_Matrix.Transpose();
+	Gauss::GaussJordan(a_Matrix, false, false);
+	a_Matrix.Transpose();
+	return {std::move(a_Matrix)};
 }
 
 // https://en.wikipedia.org/wiki/Kernel_(linear_algebra)#Computation_by_Gaussian_elimination
@@ -12,8 +14,10 @@ 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);
+	a_Matrix.Transpose();
+	a_Matrix.Augment(Matrix::Identity(a_Matrix.GetRawCount()));
+	Gauss::GaussJordan(a_Matrix, false, true);
+	a_Matrix.Transpose();
 
 	// nombre de colonnes non nulles
 	std::size_t origine_colonne = Vect(a_Matrix.SubMatrix(0, 0, matrixRawCount, matrixColumnCount)).GetCardinal();

src/Vect.cpp

@@ -2,8 +2,6 @@
 
 #include "Gauss.h"
 #include "Solver.h"
-#include <cassert>
-#include <iostream>
 
 static bool IsColumnNull(Matrix& mat, std::size_t column) {
 	for (std::size_t i = 0; i < mat.GetRawCount(); i++) {
@@ -14,7 +12,7 @@ static bool IsColumnNull(Matrix& mat, std::size_t column) {
 	return true;
 }
 
-Vect::Vect(const Matrix& a_Matrix) : m_Data(a_Matrix) {
+Vect::Vect(Matrix&& a_Matrix) : m_Data(std::move(a_Matrix)) {
 	m_Data.Transpose();
 	Gauss::GaussJordan(m_Data, false, false);
 	m_Data.Transpose();

test/test_solver.cpp

@@ -34,8 +34,8 @@ void TestKernelImage() {
 		Matrix mat, imageMat, noyauMat;
 		in >> mat >> imageMat >> noyauMat;
 
-		Vect image {imageMat};
-		Vect noyau {noyauMat};
+		Vect image {std::move(imageMat)};
+		Vect noyau {std::move(noyauMat)};
 
 		Solver solver;