diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 2a04ec4a9..8f3b7dfc1 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -26,8 +26,6 @@ #ifndef EIGEN_PARTIALLU_H #define EIGEN_PARTIALLU_H -template struct ei_partialpivlu_solve_impl; - /** \ingroup LU_Module * * \class PartialPivLU @@ -59,10 +57,11 @@ template struct ei_partialpivlu_solve_impl; * * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU */ -template class PartialPivLU +template class PartialPivLU { public: + typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef Matrix IntRowVectorType; @@ -134,11 +133,11 @@ template class PartialPivLU * \sa TriangularView::solve(), inverse(), computeInverse() */ template - inline const ei_partialpivlu_solve_impl + inline const ei_solve_return_value solve(const MatrixBase& b) const { ei_assert(m_isInitialized && "PartialPivLU is not initialized."); - return ei_partialpivlu_solve_impl(*this, b.derived()); + return ei_solve_return_value(*this, b.derived()); } /** \returns the inverse of the matrix of which *this is the LU decomposition. @@ -148,10 +147,10 @@ template class PartialPivLU * * \sa MatrixBase::inverse(), LU::inverse() */ - inline const ei_partialpivlu_solve_impl > inverse() const + inline const ei_solve_return_value > inverse() const { ei_assert(m_isInitialized && "PartialPivLU is not initialized."); - return ei_partialpivlu_solve_impl > + return ei_solve_return_value > (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); } @@ -170,6 +169,9 @@ template class PartialPivLU */ typename ei_traits::Scalar determinant() const; + inline int rows() const { return m_lu.rows(); } + inline int cols() const { return m_lu.cols(); } + protected: MatrixType m_lu; IntColVectorType m_p; @@ -407,33 +409,11 @@ typename ei_traits::Scalar PartialPivLU::determinant() c /***** Implementation of solve() *****************************************************/ -template -struct ei_traits > +template +struct ei_solve_impl, Rhs, Dest> + : ei_solve_return_value, Rhs> { - typedef Matrix ReturnMatrixType; -}; - -template -struct ei_partialpivlu_solve_impl : public ReturnByValue > -{ - typedef typename ei_cleantype::type RhsNested; - typedef PartialPivLU LUType; - const LUType& m_lu; - const typename Rhs::Nested m_rhs; - - ei_partialpivlu_solve_impl(const LUType& lu, const Rhs& rhs) - : m_lu(lu), m_rhs(rhs) - {} - - inline int rows() const { return m_lu.matrixLU().cols(); } - inline int cols() const { return m_rhs.cols(); } - - template void evalTo(Dest& dst) const + void evalTo(Dest& dst) const { /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. * So we proceed as follows: @@ -442,19 +422,22 @@ struct ei_partialpivlu_solve_impl : public ReturnByValue& dec = this->m_dec; + const Rhs& rhs = this->m_rhs; + + const int size = dec.matrixLU().rows(); + ei_assert(rhs.rows() == size); - dst.resize(size, m_rhs.cols()); + dst.resize(size, rhs.cols()); // Step 1 - for(int i = 0; i < size; ++i) dst.row(m_lu.permutationP().coeff(i)) = m_rhs.row(i); + for(int i = 0; i < size; ++i) dst.row(dec.permutationP().coeff(i)) = rhs.row(i); // Step 2 - m_lu.matrixLU().template triangularView().solveInPlace(dst); + dec.matrixLU().template triangularView().solveInPlace(dst); // Step 3 - m_lu.matrixLU().template triangularView().solveInPlace(dst); + dec.matrixLU().template triangularView().solveInPlace(dst); } };