big internal rework

src/Gauss.cpp

@@ -4,72 +4,66 @@
 
 namespace Gauss {
 
-static void GaussNonJordan(Matrix& mat, bool reduite) {
-	int r = -1;
-	for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {
-		std::size_t indice_ligne_maximum = r + 1;
-
-		// Recherche maximum
-		for (std::size_t i = r + 1; i < mat.GetRawCount(); i++) {
-			if (std::abs(mat.at(i, j)) > std::abs(mat.at(indice_ligne_maximum, j)))
-				indice_ligne_maximum = i;
-		}
-
-		// Si A[k,j]≠0 alors    (A[k,j] désigne la valeur de la ligne k et de la colonne j)
-		if (!IsEqualZero(mat.at(indice_ligne_maximum, j))) {
-			r++;
-
-			// Si k≠r alors
-			if (indice_ligne_maximum != r) {
-				// Échanger les lignes k et r  (On place la ligne du pivot en position r)
-				for (std::size_t k = 0; k < mat.GetColumnCount(); k++) {
-					std::swap(mat.at(indice_ligne_maximum, k), mat.at(r, k));
-				}
-			}
-
-			// Pour i de 1 jusqu'à n               (On simplifie les autres lignes)
-			for (std::size_t i = (reduite ? 0 : j); i < mat.GetRawCount(); i++) {
-				// Si i≠r alors
-				if (i != r) {
-					// Soustraire à la ligne i la ligne r multipliée par A[i,j] (de façon à
-					// annuler A[i,j])
-					for (int k = mat.GetColumnCount() - 1; k >= 0; k--) {
-						long double pivot = mat.at(r, j);
-						long double anul = mat.at(i, j);
-						mat.at(i, k) = mat.at(i, k) * pivot - mat.at(r, k) * anul;
-					}
-				}
-			}
-		}
+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 GaussJordan(Matrix& mat, bool reduite) {
-	GaussNonJordan(mat, reduite);
-	for (std::size_t i = 0; i < mat.GetRawCount(); i++) {
-		int k = -1;
-		for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {
-			if (!IsEqualZero(mat.at(i, j))) {
-				k = j;
-				break;
-			}
-		}
-		// ligne de 0
-		if (k == -1)
-			break;
-		// on divise la ligne par (i, k)
-		long double annul = mat.at(i, k);
-		for (int j = 0; j < mat.GetColumnCount(); j++) {
-			mat.at(i, j) /= annul;
+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 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))) {
+			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++) {
+		mat.at(line, j) = mat.at(line, j) * pivot - mat.at(pivot_line, j) * anul;
+	}
 }
 
 void GaussJordan(Matrix& mat, bool reduite, bool normalise) {
-	if (normalise)
-		GaussJordan(mat, reduite);
-	else
-		GaussNonJordan(mat, reduite);
+	int indice_ligne_pivot = -1;
+
+	for (std::size_t j = 0; j < mat.GetColumnCount(); j++) {
+
+		int indice_ligne_pivot_trouve = FirstNotNullElementIndexOnColumn(mat, 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(mat, indice_ligne_pivot_trouve, indice_ligne_pivot);
+		}
+
+		Matrix::Element pivot = mat.at(indice_ligne_pivot, j);
+
+		if (normalise) {
+			DivideLine(mat, indice_ligne_pivot, pivot);
+		}
+
+		// On simplifie les autres lignes
+		for (std::size_t i = (reduite ? 0 : j); i < mat.GetRawCount(); i++) {
+			// Pour les lignes autre que la ligne pivot
+			if (i != static_cast<std::size_t>(indice_ligne_pivot)) {
+				SimplifyLine(mat, i, indice_ligne_pivot, j);
+			}
+		}
+	}
 }
 
 }  // namespace Gauss