From 5391b7b76a9a2862a9f25c3384d4cced61ba7112 Mon Sep 17 00:00:00 2001 From: Persson-dev Date: Thu, 9 May 2024 15:56:59 +0200 Subject: [PATCH] more optimizations --- include/Gauss.h | 8 ---- include/Matrix.h | 6 +++ include/Vect.h | 2 +- src/Gauss.cpp | 102 ++++++++++++++++++------------------------- src/Matrix.cpp | 12 +++++ src/Solver.cpp | 12 +++-- src/Vect.cpp | 4 +- test/test_solver.cpp | 4 +- 8 files changed, 72 insertions(+), 78 deletions(-) diff --git a/include/Gauss.h b/include/Gauss.h index 0a50ec5..a21bfb3 100644 --- a/include/Gauss.h +++ b/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 \ No newline at end of file diff --git a/include/Matrix.h b/include/Matrix.h index 816b21b..494e003 100644 --- a/include/Matrix.h +++ b/include/Matrix.h @@ -17,6 +17,7 @@ class Matrix { public: typedef long double Element; + typedef std::vector::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&& a_Elements); + + iterator begin(); + iterator end(); + + iterator GetLineIterator(std::size_t a_Raw); }; template diff --git a/include/Vect.h b/include/Vect.h index 437551a..ccaa425 100644 --- a/include/Vect.h +++ b/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 diff --git a/src/Gauss.cpp b/src/Gauss.cpp index dff05c6..1bdebee 100644 --- a/src/Gauss.cpp +++ b/src/Gauss.cpp @@ -2,30 +2,20 @@ #include "Matrix.h" +#include +#include +#include + 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(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(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(indice_ligne_pivot)) { SimplifyLine(a_Matrix, i, indice_ligne_pivot, j); } - } + }); } } -void GaussJordanColumn(Matrix& a_Matrix, bool a_Reduite, bool a_Normalise) { - int indice_colonne_pivot = -1; +static void GaussJordanTriangular(Matrix& a_Matrix, bool a_Normalise) { + int indice_ligne_pivot = -1; - for (std::size_t j = 0; j < a_Matrix.GetRawCount(); j++) { + for (std::size_t j = 0; j < a_Matrix.GetColumnCount(); j++) { - int indice_colonne_pivot_trouve = FirstNotNullElementIndexOnLine(a_Matrix, j, indice_colonne_pivot + 1); + int indice_ligne_pivot_trouve = FirstNotNullElementIndexOnColumn(a_Matrix, j, indice_ligne_pivot + 1); - if (indice_colonne_pivot_trouve < 0) // ligne de 0 - continue; // on regarde la prochaine ligne + if (indice_ligne_pivot_trouve < 0) // colonne de 0 + continue; // on regarde la prochaine colonne - indice_colonne_pivot++; + indice_ligne_pivot++; - if (indice_colonne_pivot_trouve != indice_colonne_pivot) { - SwapColumns(a_Matrix, indice_colonne_pivot_trouve, indice_colonne_pivot); + if (indice_ligne_pivot_trouve != indice_ligne_pivot) { + SwapLines(a_Matrix, indice_ligne_pivot_trouve, indice_ligne_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(indice_colonne_pivot)) { - SimplifyColumn(a_Matrix, i, j, indice_colonne_pivot); - } - } + auto range = std::views::iota(static_cast(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 \ No newline at end of file diff --git a/src/Matrix.cpp b/src/Matrix.cpp index 5fbd3e1..de8b9f2 100644 --- a/src/Matrix.cpp +++ b/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(); +} \ No newline at end of file diff --git a/src/Solver.cpp b/src/Solver.cpp index f4c35f1..50d0a02 100644 --- a/src/Solver.cpp +++ b/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(); diff --git a/src/Vect.cpp b/src/Vect.cpp index c1b29d8..4747ae3 100644 --- a/src/Vect.cpp +++ b/src/Vect.cpp @@ -2,8 +2,6 @@ #include "Gauss.h" #include "Solver.h" -#include -#include 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(); diff --git a/test/test_solver.cpp b/test/test_solver.cpp index 83a93bb..83f1a07 100644 --- a/test/test_solver.cpp +++ b/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;