From 79a0feee683f7eef84134c85416e52519ffaf350 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Mon, 4 Aug 2008 23:34:21 +0000
Subject: [PATCH] big performance improvement in inverse and LU

---
 Eigen/src/Core/Matrix.h |  1 -
 Eigen/src/LU/Inverse.h  | 17 ++++++++---------
 Eigen/src/LU/LU.h       |  4 ++--
 3 files changed, 10 insertions(+), 12 deletions(-)

diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index 85a984920..a1c313108 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -179,7 +179,6 @@ class Matrix : public MatrixBase<Matrix<_Scalar, _Rows, _Cols, _MaxRows, _MaxCol
       ei_pstoret<Scalar, PacketScalar, StoreMode>(m_storage.data() + index, x);
     }
 
-  public:
     /** \returns a const pointer to the data array of this matrix */
     inline const Scalar *data() const
     { return m_storage.data(); }
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index 68cac5985..1a1fc2334 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -50,10 +50,10 @@ struct ei_compute_inverse_in_general_case<MatrixType, RowMajor>
       matrix.row(k).swap(matrix.row(k+rowOfBiggest));
 
       const Scalar d = matrix(k,k);
-      inverse.block(k+1, 0, size-k-1, size)
-        -= matrix.col(k).end(size-k-1) * (inverse.row(k) / d);
-      matrix.corner(BottomRight, size-k-1, size-k)
-        -= matrix.col(k).end(size-k-1) * (matrix.row(k).end(size-k) / d);
+      for(int row = k + 1; row < size; row++)
+        inverse.row(row) -= inverse.row(k) * (matrix.coeff(row,k)/d);
+      for(int row = k + 1; row < size; row++)
+        matrix.row(row).end(size-k) -= (matrix.row(k).end(size-k)/d) * matrix.coeff(row,k);
     }
 
     for(int k = 0; k < size-1; k++)
@@ -87,12 +87,11 @@ struct ei_compute_inverse_in_general_case<MatrixType, ColMajor>
       matrix.row(k).end(size-k).cwise().abs().maxCoeff(&colOfBiggest);
       inverse.col(k).swap(inverse.col(k+colOfBiggest));
       matrix.col(k).swap(matrix.col(k+colOfBiggest));
-
       const Scalar d = matrix(k,k);
-      inverse.block(0, k+1, size, size-k-1)
-        -= (inverse.col(k) / d) * matrix.row(k).end(size-k-1);
-      matrix.corner(BottomRight, size-k, size-k-1)
-        -= (matrix.col(k).end(size-k) / d) * matrix.row(k).end(size-k-1);
+      for(int col = k + 1; col < size; col++)
+        inverse.col(col) -= inverse.col(k) * (matrix.coeff(k,col)/d);
+      for(int col = k + 1; col < size; col++)
+        matrix.col(col).end(size-k) -= (matrix.col(k).end(size-k)/d) * matrix.coeff(k,col);
     }
 
     for(int k = 0; k < size-1; k++)
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index aa701f2ce..4da026dda 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -143,8 +143,8 @@ LU<MatrixType>::LU(const MatrixType& matrix)
     if(k<rows-1)
       m_lu.col(k).end(rows-k-1) /= lu_k_k;
     if(k<size-1)
-      m_lu.corner(BottomRight, rows-k-1, cols-k-1)
-          -= m_lu.col(k).end(rows-k-1) * m_lu.row(k).end(cols-k-1);
+      for( int col = k + 1; col < cols; col++ )
+        m_lu.col(col).end(rows-k-1) -= m_lu.col(k).end(rows-k-1) * m_lu.coeff(k,col);
   }
 
   for(int k = 0; k < matrix.rows(); k++) m_p.coeffRef(k) = k;