Merge pull request #156 from SimonRohou/codac2_dev

[mat] corrected bug in Gauss Jordan decomposition (thanks Damien Massé!)
codac-team · Dec 9, 2024 · 0afe027 · 0afe027
2 parents 639b848 + 6875f2f
commit 0afe027
Show file tree

Hide file tree

Showing 3 changed files with 13 additions and 10 deletions.
diff --git a/src/core/matrices/codac2_GaussJordan.cpp b/src/core/matrices/codac2_GaussJordan.cpp
@@ -2,7 +2,7 @@
  *  codac2_GaussJordan.cpp
  * ----------------------------------------------------------------------------
  *  \date       2024
- *  \author     Simon Rohou
+ *  \author     Luc Jaulin, Simon Rohou, Damien Massé
  *  \copyright  Copyright 2024 Codac Team
  *  \license    GNU Lesser General Public License (LGPL)
  */
@@ -40,17 +40,18 @@ namespace codac2
 
   Matrix gauss_jordan(const Matrix& A)
   {
-    Index n = A.rows();//, m = A.cols();
+    Index n = A.rows(), m = A.cols();
     Eigen::FullPivLU<Matrix> lu(A);
 
     Matrix L = Matrix::Identity(n,n);
-    //if(std::pow(L.determinant(),2) < 1e-5)
-    //{
-    //  cout << "[Matrix gauss_jordan(const Matrix& A)] -> eye Matrix" << endl;
-    //  return Matrix::eye(n,n);
-    //}
+    if(std::pow(L.determinant(),2) < 1e-5)
+    {
+      cout << "[Matrix gauss_jordan(const Matrix& A)] -> eye Matrix" << endl;
+      return Matrix::eye(n,n);
+    }
 
-    L.block(0,0,n,n).triangularView<Eigen::StrictlyLower>() = lu.matrixLU();
+    L.block(0,0,n,std::min(m,n)).triangularView<Eigen::StrictlyLower>() = 
+        lu.matrixLU().block(0,0,n,std::min(m,n));
 
     Matrix P = lu.permutationP();
     Matrix U = lu.matrixLU().triangularView<Eigen::Upper>();

diff --git a/src/core/matrices/codac2_GaussJordan.h b/src/core/matrices/codac2_GaussJordan.h
@@ -2,7 +2,7 @@
  *  \file codac2_GaussJordan.h
  * ----------------------------------------------------------------------------
  *  \date       2024
- *  \author     Simon Rohou
+ *  \author     Luc Jaulin, Simon Rohou, Damien Massé
  *  \copyright  Copyright 2024 Codac Team
  *  \license    GNU Lesser General Public License (LGPL)
  */
@@ -16,4 +16,6 @@ namespace codac2
 {
   // Gauss Jordan band diagonalization preconditioning
   Matrix gauss_jordan(const Matrix& A);
+
+  // From https://www.ensta-bretagne.fr/jaulin/centered.html
 }
diff --git a/src/core/matrices/codac2_matrices.h b/src/core/matrices/codac2_matrices.h
@@ -41,7 +41,7 @@ namespace Eigen
  * (this is a standard C++ macro which disables all asserts).
  * https://eigen.tuxfamily.org/dox/TopicPreprocessorDirectives.html
  */
-#define EIGEN_NO_DEBUG
+//#define EIGEN_NO_DEBUG // uncomment to disable Eigen's assertions
 #endif
 
 #include <Eigen/Core>