diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-12-22 22:51:08 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-12-22 22:51:08 +0100 |
commit | eaaba30cacf078335ad41314bf1e4ce993f0b3b4 (patch) | |
tree | 61156273b56dc19fc4e661a5aeb42d8e6112a307 /Eigen/src/LU/PartialPivLU.h | |
parent | e182e9c6166840b8db44e0be459a2e26d26fda0f (diff) | |
parent | eeec7d842e05394703c42e7569230835f7842089 (diff) |
merge with default branch
Diffstat (limited to 'Eigen/src/LU/PartialPivLU.h')
-rw-r--r-- | Eigen/src/LU/PartialPivLU.h | 29 |
1 files changed, 12 insertions, 17 deletions
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 37cffd7a1..deb29b5c6 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -64,10 +64,8 @@ template<typename _MatrixType> class PartialPivLU typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; - typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; - typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType; + typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> PermutationVectorType; + typedef PermutationMatrix<MatrixType::RowsAtCompileTime> PermutationType; enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN( MatrixType::MaxColsAtCompileTime, @@ -105,11 +103,9 @@ template<typename _MatrixType> class PartialPivLU return m_lu; } - /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, - * representing the P permutation i.e. the permutation of the rows. For its precise meaning, - * see the examples given in the documentation of class FullPivLU. + /** \returns the permutation matrix P. */ - inline const IntColVectorType& permutationP() const + inline const PermutationType& permutationP() const { ei_assert(m_isInitialized && "PartialPivLU is not initialized."); return m_p; @@ -147,11 +143,11 @@ template<typename _MatrixType> class PartialPivLU * * \sa MatrixBase::inverse(), LU::inverse() */ - inline const ei_solve_retval<PartialPivLU,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const + inline const ei_solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const { ei_assert(m_isInitialized && "PartialPivLU is not initialized."); - return ei_solve_retval<PartialPivLU,NestByValue<typename MatrixType::IdentityReturnType> > - (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); + return ei_solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> + (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols())); } /** \returns the determinant of the matrix of which @@ -174,7 +170,7 @@ template<typename _MatrixType> class PartialPivLU protected: MatrixType m_lu; - IntColVectorType m_p; + PermutationType m_p; int m_det_p; bool m_isInitialized; }; @@ -379,20 +375,19 @@ template<typename MatrixType> PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix) { m_lu = matrix; - m_p.resize(matrix.rows()); ei_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); const int size = matrix.rows(); - IntColVectorType rows_transpositions(size); + PermutationVectorType rows_transpositions(size); int nb_transpositions; ei_partial_lu_inplace(m_lu, rows_transpositions, nb_transpositions); m_det_p = (nb_transpositions%2) ? -1 : 1; - for(int k = 0; k < size; ++k) m_p.coeffRef(k) = k; + m_p.setIdentity(size); for(int k = size-1; k >= 0; --k) - std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k))); + m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k)); m_isInitialized = true; return *this; @@ -428,7 +423,7 @@ struct ei_solve_retval<PartialPivLU<_MatrixType>, Rhs> dst.resize(size, rhs().cols()); // Step 1 - for(int i = 0; i < size; ++i) dst.row(dec().permutationP().coeff(i)) = rhs().row(i); + dst = dec().permutationP() * rhs(); // Step 2 dec().matrixLU().template triangularView<UnitLowerTriangular>().solveInPlace(dst); |