From 09ef7db9d92f6f96ca90655ae3ceecbefa4421fd Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Tue, 5 Aug 2008 15:43:11 +0000 Subject: [PATCH] Add partial pivoting runtime option to LU. Note: in fact, inverse() always uses partial pivoting because the algo currently used doesn't make sense with complete pivoting. No num stability issue so far even with size 200x200. If there is any problem we can of course reimplement inverse on top of LU. --- Eigen/src/Core/MatrixBase.h | 2 +- Eigen/src/Core/util/Constants.h | 5 ++ Eigen/src/LU/LU.h | 98 +++++++++++++++++++++++---------- 3 files changed, 74 insertions(+), 31 deletions(-) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index c2b12a0a3..4b2c26341 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -526,7 +526,7 @@ template class MatrixBase /////////// LU module /////////// - const LU lu() const; + const LU lu(int pivoting) const; const EvalType inverse() const; void computeInverse(EvalType *result) const; Scalar determinant() const; diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 24c653e2e..dd854e9ac 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -209,6 +209,11 @@ enum { HasDirectAccess = DirectAccessBit }; +enum { + PartialPivoting, + CompletePivoting +}; + const int FullyCoherentAccessPattern = 0x1; const int InnerCoherentAccessPattern = 0x2 | FullyCoherentAccessPattern; const int OuterCoherentAccessPattern = 0x4 | InnerCoherentAccessPattern; diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index 4da026dda..a707ce8d7 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -54,29 +54,29 @@ template class LU typedef Matrix IntRowVectorType; typedef Matrix IntColVectorType; - LU(const MatrixType& matrix); + LU(const MatrixType& matrix, int pivoting = CompletePivoting); - const MatrixType& matrixLU() const + inline const MatrixType& matrixLU() const { return m_lu; } - const Part matrixL() const + inline const Part matrixL() const { return m_lu; } - const Part matrixU() const + inline const Part matrixU() const { return m_lu; } - const IntColVectorType& permutationP() const + inline const IntColVectorType& permutationP() const { return m_p; } - const IntRowVectorType& permutationQ() const + inline const IntRowVectorType& permutationQ() const { return m_q; } @@ -92,29 +92,47 @@ template class LU * as the LU decomposition has already been computed. * * Warning: a determinant can be very big or small, so for matrices - * of large dimension (like a 50-by-50 matrix) there can be a risk of + * of large enough dimension (like a 50-by-50 matrix) there is a risk of * overflow/underflow. */ typename ei_traits::Scalar determinant() const; + inline int rank() const + { + return m_rank; + } + + inline int dimensionOfKernel() const + { + return m_lu.cols() - m_rank; + } + + inline bool isInvertible() const + { + return m_rank == m_lu.cols(); + } + protected: MatrixType m_lu; IntColVectorType m_p; IntRowVectorType m_q; int m_det_pq; Scalar m_biggest_eigenvalue_of_u; - int m_dimension_of_kernel; + int m_rank; + int m_pivoting; }; template -LU::LU(const MatrixType& matrix) +LU::LU(const MatrixType& matrix, int pivoting) : m_lu(matrix), m_p(matrix.rows()), - m_q(matrix.cols()) + m_q(matrix.cols()), + m_pivoting(pivoting) { const int size = matrix.diagonal().size(); const int rows = matrix.rows(); const int cols = matrix.cols(); + ei_assert(pivoting == PartialPivoting || pivoting == CompletePivoting); IntColVectorType rows_transpositions(matrix.rows()); IntRowVectorType cols_transpositions(matrix.cols()); @@ -123,20 +141,36 @@ LU::LU(const MatrixType& matrix) for(int k = 0; k < size; k++) { int row_of_biggest, col_of_biggest; - const Scalar biggest = m_lu.corner(Eigen::BottomRight, rows-k, cols-k) - .cwise().abs() - .maxCoeff(&row_of_biggest, &col_of_biggest); - row_of_biggest += k; - col_of_biggest += k; - rows_transpositions.coeffRef(k) = row_of_biggest; - cols_transpositions.coeffRef(k) = col_of_biggest; - if(k != row_of_biggest) { - m_lu.row(k).swap(m_lu.row(row_of_biggest)); - number_of_transpositions++; + Scalar biggest; + if(m_pivoting == CompletePivoting) + { + biggest = m_lu.corner(Eigen::BottomRight, rows-k, cols-k) + .cwise().abs() + .maxCoeff(&row_of_biggest, &col_of_biggest); + row_of_biggest += k; + col_of_biggest += k; + rows_transpositions.coeffRef(k) = row_of_biggest; + cols_transpositions.coeffRef(k) = col_of_biggest; + if(k != row_of_biggest) { + m_lu.row(k).swap(m_lu.row(row_of_biggest)); + number_of_transpositions++; + } + if(k != col_of_biggest) { + m_lu.col(k).swap(m_lu.col(col_of_biggest)); + number_of_transpositions++; + } } - if(k != col_of_biggest) { - m_lu.col(k).swap(m_lu.col(col_of_biggest)); - number_of_transpositions++; + else // partial pivoting + { + biggest = m_lu.col(k).end(rows-k) + .cwise().abs() + .maxCoeff(&row_of_biggest); + row_of_biggest += k; + rows_transpositions.coeffRef(k) = row_of_biggest; + if(k != row_of_biggest) { + m_lu.row(k).swap(m_lu.row(row_of_biggest)); + number_of_transpositions++; + } } const Scalar lu_k_k = m_lu.coeff(k,k); if(ei_isMuchSmallerThan(lu_k_k, biggest)) continue; @@ -151,9 +185,12 @@ LU::LU(const MatrixType& matrix) for(int k = size-1; k >= 0; k--) std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k))); - for(int k = 0; k < matrix.cols(); k++) m_q.coeffRef(k) = k; - for(int k = 0; k < size; k++) - std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k))); + if(pivoting == CompletePivoting) + { + for(int k = 0; k < matrix.cols(); k++) m_q.coeffRef(k) = k; + for(int k = 0; k < size; k++) + std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k))); + } m_det_pq = (number_of_transpositions%2) ? -1 : 1; @@ -161,14 +198,15 @@ LU::LU(const MatrixType& matrix) m_lu.diagonal().cwise().abs().maxCoeff(&index_of_biggest); m_biggest_eigenvalue_of_u = m_lu.diagonal().coeff(index_of_biggest); - m_dimension_of_kernel = 0; + m_rank = 0; for(int k = 0; k < size; k++) - m_dimension_of_kernel += ei_isMuchSmallerThan(m_lu.diagonal().coeff(k), m_biggest_eigenvalue_of_u); + m_rank += !ei_isMuchSmallerThan(m_lu.diagonal().coeff(k), m_biggest_eigenvalue_of_u); } template typename ei_traits::Scalar LU::determinant() const { + if(!isInvertible()) return Scalar(0); Scalar res = m_det_pq; for(int k = 0; k < m_lu.diagonal().size(); k++) res *= m_lu.diagonal().coeff(k); return res; @@ -180,9 +218,9 @@ typename ei_traits::Scalar LU::determinant() const */ template const LU::EvalType> -MatrixBase::lu() const +MatrixBase::lu(int pivoting = CompletePivoting) const { - return eval(); + return LU::EvalType>(eval(), pivoting); } #endif // EIGEN_LU_H