diff options
Diffstat (limited to 'Eigen/src/LU/PartialPivLU.h')
-rw-r--r-- | Eigen/src/LU/PartialPivLU.h | 85 |
1 files changed, 70 insertions, 15 deletions
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 1d389ecac..ac53f7ab5 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -13,6 +13,19 @@ namespace Eigen { +namespace internal { +template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> > + : traits<_MatrixType> +{ + typedef traits<_MatrixType> BaseTraits; + enum { + Flags = BaseTraits::Flags & RowMajorBit, + CoeffReadCost = Dynamic + }; +}; + +} // end namespace internal + /** \ingroup LU_Module * * \class PartialPivLU @@ -62,6 +75,7 @@ template<typename _MatrixType> class PartialPivLU typedef typename MatrixType::Index Index; typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType; typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; + typedef typename MatrixType::PlainObject PlainObject; /** @@ -128,6 +142,15 @@ template<typename _MatrixType> class PartialPivLU * * \sa TriangularView::solve(), inverse(), computeInverse() */ +#ifdef EIGEN_TEST_EVALUATORS + template<typename Rhs> + inline const Solve<PartialPivLU, Rhs> + solve(const MatrixBase<Rhs>& b) const + { + eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); + return Solve<PartialPivLU, Rhs>(*this, b.derived()); + } +#else template<typename Rhs> inline const internal::solve_retval<PartialPivLU, Rhs> solve(const MatrixBase<Rhs>& b) const @@ -135,6 +158,7 @@ template<typename _MatrixType> class PartialPivLU eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); return internal::solve_retval<PartialPivLU, Rhs>(*this, b.derived()); } +#endif /** \returns the inverse of the matrix of which *this is the LU decomposition. * @@ -143,12 +167,20 @@ template<typename _MatrixType> class PartialPivLU * * \sa MatrixBase::inverse(), LU::inverse() */ +#ifdef EIGEN_TEST_EVALUATORS + inline const Inverse<PartialPivLU> inverse() const + { + eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); + return Inverse<PartialPivLU>(*this); + } +#else inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const { eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols())); } +#endif /** \returns the determinant of the matrix of which * *this is the LU decomposition. It has only linear complexity @@ -169,6 +201,30 @@ template<typename _MatrixType> class PartialPivLU inline Index rows() const { return m_lu.rows(); } inline Index cols() const { return m_lu.cols(); } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template<typename RhsType, typename DstType> + EIGEN_DEVICE_FUNC + void _solve_impl(const RhsType &rhs, DstType &dst) const { + /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. + * So we proceed as follows: + * Step 1: compute c = Pb. + * Step 2: replace c by the solution x to Lx = c. + * Step 3: replace c by the solution x to Ux = c. + */ + + eigen_assert(rhs.rows() == m_lu.rows()); + + // Step 1 + dst = permutationP() * rhs; + + // Step 2 + m_lu.template triangularView<UnitLower>().solveInPlace(dst); + + // Step 3 + m_lu.template triangularView<Upper>().solveInPlace(dst); + } + #endif protected: MatrixType m_lu; @@ -434,6 +490,7 @@ MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const namespace internal { +#ifndef EIGEN_TEST_EVALUATORS template<typename _MatrixType, typename Rhs> struct solve_retval<PartialPivLU<_MatrixType>, Rhs> : solve_retval_base<PartialPivLU<_MatrixType>, Rhs> @@ -442,23 +499,21 @@ struct solve_retval<PartialPivLU<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. - * So we proceed as follows: - * Step 1: compute c = Pb. - * Step 2: replace c by the solution x to Lx = c. - * Step 3: replace c by the solution x to Ux = c. - */ - - eigen_assert(rhs().rows() == dec().matrixLU().rows()); - - // Step 1 - dst = dec().permutationP() * rhs(); + dec()._solve_impl(rhs(), dst); + } +}; +#endif - // Step 2 - dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst); +/***** Implementation of inverse() *****************************************************/ - // Step 3 - dec().matrixLU().template triangularView<Upper>().solveInPlace(dst); +template<typename DstXprType, typename MatrixType, typename Scalar> +struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar> +{ + typedef PartialPivLU<MatrixType> LuType; + typedef Inverse<LuType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &) + { + dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); } }; |