From becbeda50ac17288dba0a93c6adc67b663d32a7a Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sat, 9 Aug 2008 15:50:07 +0000 Subject: [PATCH] * reimplement the general case of inverse() on top of LU. Advantages: - removes much code - 2.5x faster (even though LU uses complete pivoting contrary to what inverse used to do!) - there _were_ numeric stability problems with the partial pivoting approach of inverse(), with 200x200 matrices inversion failed almost half of the times (overflow). Now these problems are solved thanks to complete pivoting. * remove some useless stuff in LU --- Eigen/src/LU/Inverse.h | 90 ++---------------------------------------- Eigen/src/LU/LU.h | 24 +++-------- 2 files changed, 9 insertions(+), 105 deletions(-) diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index dd34d65e2..6f6256d92 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -25,90 +25,8 @@ #ifndef EIGEN_INVERSE_H #define EIGEN_INVERSE_H -/*************************************************************************** -*** Part 1 : implementation in the general case, by Gaussian elimination *** -***************************************************************************/ - -template -struct ei_compute_inverse_in_general_case; - -template -struct ei_compute_inverse_in_general_case -{ - static void run(const MatrixType& _matrix, MatrixType *result) - { - typedef typename MatrixType::Scalar Scalar; - MatrixType matrix(_matrix); - MatrixType &inverse = *result; - inverse.setIdentity(); - const int size = matrix.rows(); - for(int k = 0; k < size-1; k++) - { - int rowOfBiggest; - matrix.col(k).end(size-k).cwise().abs().maxCoeff(&rowOfBiggest); - inverse.row(k).swap(inverse.row(k+rowOfBiggest)); - matrix.row(k).swap(matrix.row(k+rowOfBiggest)); - - const Scalar d = matrix(k,k); - 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++) - { - const Scalar d = static_cast(1)/matrix(k,k); - matrix.row(k).end(size-k) *= d; - inverse.row(k) *= d; - } - inverse.row(size-1) /= matrix(size-1,size-1); - - for(int k = size-1; k >= 1; k--) - for(int row = 0; row < k; row++) - inverse.row(row) -= inverse.row(k) * matrix.coeff(row,k); - } -}; - -template -struct ei_compute_inverse_in_general_case -{ - static void run(const MatrixType& _matrix, MatrixType *result) - { - typedef typename MatrixType::Scalar Scalar; - MatrixType matrix(_matrix); - MatrixType& inverse = *result; - inverse.setIdentity(); - const int size = matrix.rows(); - for(int k = 0; k < size-1; k++) - { - int colOfBiggest; - 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); - 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++) - { - const Scalar d = static_cast(1)/matrix(k,k); - matrix.col(k).end(size-k) *= d; - inverse.col(k) *= d; - } - inverse.col(size-1) /= matrix(size-1,size-1); - - for(int k = size-1; k >= 1; k--) - for(int col = 0; col < k; col++) - inverse.col(col) -= inverse.col(k) * matrix.coeff(k,col); - } -}; - /******************************************************************** -*** Part 2 : optimized implementations for fixed-size 2,3,4 cases *** +*** Part 1 : optimized implementations for fixed-size 2,3,4 cases *** ********************************************************************/ template @@ -234,7 +152,7 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu } /*********************************************** -*** Part 3 : selector and MatrixBase methods *** +*** Part 2 : selector and MatrixBase methods *** ***********************************************/ template @@ -242,8 +160,8 @@ struct ei_compute_inverse { static inline void run(const MatrixType& matrix, MatrixType* result) { - ei_compute_inverse_in_general_case - ::run(matrix, result); + LU lu(matrix); + lu.solve(MatrixType::Identity(matrix.rows(), matrix.cols()), result); } }; diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index 6a55921fe..e603d09dc 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -58,16 +58,7 @@ template class LU enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN( MatrixType::MaxColsAtCompileTime, - MatrixType::MaxRowsAtCompileTime), - SmallDimAtCompileTime = EIGEN_ENUM_MIN( - MatrixType::ColsAtCompileTime, - MatrixType::RowsAtCompileTime), - MaxBigDimAtCompileTime = EIGEN_ENUM_MAX( - MatrixType::MaxColsAtCompileTime, - MatrixType::MaxRowsAtCompileTime), - BigDimAtCompileTime = EIGEN_ENUM_MAX( - MatrixType::ColsAtCompileTime, - MatrixType::RowsAtCompileTime) + MatrixType::MaxRowsAtCompileTime) }; LU(const MatrixType& matrix); @@ -107,12 +98,10 @@ template class LU MatrixType::MaxColsAtCompileTime, LU::MaxSmallDimAtCompileTime> kernel() const; - template + template bool solve( const MatrixBase& b, - Matrix *result + ResultType *result ) const; /** @@ -281,12 +270,10 @@ LU::kernel() const } template -template +template bool LU::solve( const MatrixBase& b, - Matrix *result + ResultType *result ) const { /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}. @@ -298,7 +285,6 @@ bool LU::solve( */ ei_assert(b.rows() == m_lu.rows()); - const int bigdim = std::max(m_lu.rows(), m_lu.cols()); const int smalldim = std::min(m_lu.rows(), m_lu.cols()); typename OtherDerived::Eval c(b.rows(), b.cols());