diff --git a/Eigen/src/Array/CwiseOperators.h b/Eigen/src/Array/CwiseOperators.h index 7a8b9935b..1cd1866e7 100644 --- a/Eigen/src/Array/CwiseOperators.h +++ b/Eigen/src/Array/CwiseOperators.h @@ -43,38 +43,6 @@ Cwise::sqrt() const return _expression(); } -/** \array_module - * - * \returns an expression of the coefficient-wise exponential of *this. - * - * Example: \include Cwise_exp.cpp - * Output: \verbinclude Cwise_exp.out - * - * \sa pow(), log(), sin(), cos() - */ -template -inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_exp_op) -Cwise::exp() const -{ - return _expression(); -} - -/** \array_module - * - * \returns an expression of the coefficient-wise logarithm of *this. - * - * Example: \include Cwise_log.cpp - * Output: \verbinclude Cwise_log.out - * - * \sa exp() - */ -template -inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_log_op) -Cwise::log() const -{ - return _expression(); -} - /** \array_module * * \returns an expression of the coefficient-wise cosine of *this. diff --git a/Eigen/src/Array/Functors.h b/Eigen/src/Array/Functors.h index 53a9019a2..fd259f7bc 100644 --- a/Eigen/src/Array/Functors.h +++ b/Eigen/src/Array/Functors.h @@ -69,40 +69,6 @@ struct ei_functor_traits > }; }; -/** \internal - * - * \array_module - * - * \brief Template functor to compute the exponential of a scalar - * - * \sa class CwiseUnaryOp, Cwise::exp() - */ -template struct ei_scalar_exp_op EIGEN_EMPTY_STRUCT { - inline const Scalar operator() (const Scalar& a) const { return ei_exp(a); } - typedef typename ei_packet_traits::type Packet; - inline Packet packetOp(const Packet& a) const { return ei_pexp(a); } -}; -template -struct ei_functor_traits > -{ enum { Cost = 5 * NumTraits::MulCost, PacketAccess = ei_packet_traits::HasExp }; }; - -/** \internal - * - * \array_module - * - * \brief Template functor to compute the logarithm of a scalar - * - * \sa class CwiseUnaryOp, Cwise::log() - */ -template struct ei_scalar_log_op EIGEN_EMPTY_STRUCT { - inline const Scalar operator() (const Scalar& a) const { return ei_log(a); } - typedef typename ei_packet_traits::type Packet; - inline Packet packetOp(const Packet& a) const { return ei_plog(a); } -}; -template -struct ei_functor_traits > -{ enum { Cost = 5 * NumTraits::MulCost, PacketAccess = ei_packet_traits::HasLog }; }; - /** \internal * * \array_module diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h index 3ffb24833..d701c25d9 100644 --- a/Eigen/src/Core/CwiseUnaryOp.h +++ b/Eigen/src/Core/CwiseUnaryOp.h @@ -205,6 +205,35 @@ MatrixBase::cast() const return derived(); } +/** \returns an expression of the coefficient-wise exponential of *this. + * + * Example: \include Cwise_exp.cpp + * Output: \verbinclude Cwise_exp.out + * + * \sa pow(), log(), sin(), cos() + */ +template +inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_exp_op) +Cwise::exp() const +{ + return _expression(); +} + +/** \returns an expression of the coefficient-wise logarithm of *this. + * + * Example: \include Cwise_log.cpp + * Output: \verbinclude Cwise_log.out + * + * \sa exp() + */ +template +inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_log_op) +Cwise::log() const +{ + return _expression(); +} + + /** \relates MatrixBase */ template EIGEN_STRONG_INLINE const typename MatrixBase::ScalarMultipleReturnType diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index a4c9604df..0c68d7434 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -298,6 +298,36 @@ template struct ei_functor_traits > { enum { Cost = 0, PacketAccess = false }; }; +/** \internal + * + * \brief Template functor to compute the exponential of a scalar + * + * \sa class CwiseUnaryOp, Cwise::exp() + */ +template struct ei_scalar_exp_op EIGEN_EMPTY_STRUCT { + inline const Scalar operator() (const Scalar& a) const { return ei_exp(a); } + typedef typename ei_packet_traits::type Packet; + inline Packet packetOp(const Packet& a) const { return ei_pexp(a); } +}; +template +struct ei_functor_traits > +{ enum { Cost = 5 * NumTraits::MulCost, PacketAccess = ei_packet_traits::HasExp }; }; + +/** \internal + * + * \brief Template functor to compute the logarithm of a scalar + * + * \sa class CwiseUnaryOp, Cwise::log() + */ +template struct ei_scalar_log_op EIGEN_EMPTY_STRUCT { + inline const Scalar operator() (const Scalar& a) const { return ei_log(a); } + typedef typename ei_packet_traits::type Packet; + inline Packet packetOp(const Packet& a) const { return ei_plog(a); } +}; +template +struct ei_functor_traits > +{ enum { Cost = 5 * NumTraits::MulCost, PacketAccess = ei_packet_traits::HasLog }; }; + /** \internal * \brief Template functor to multiply a scalar by a fixed other one * diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h index 0ffcfe88c..77a0abedc 100644 --- a/Eigen/src/QR/FullPivotingHouseholderQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -37,10 +37,10 @@ * * This class performs a rank-revealing QR decomposition using Householder transformations. * - * This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal - * numerical stability. + * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal + * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivotingHouseholderQR. * - * \sa MatrixBase::householderRrqr() + * \sa MatrixBase::fullPivotingHouseholderQr() */ template class FullPivotingHouseholderQR { @@ -125,11 +125,26 @@ template class FullPivotingHouseholderQR * * \warning a determinant can be very big or small, so for matrices * of large enough dimension, there is a risk of overflow/underflow. + * One way to work around that is to use logAbsDeterminant() instead. * - * \sa MatrixBase::determinant() + * \sa logAbsDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar absDeterminant() const; + /** \returns the natural log of the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \note This method is useful to work around the risk of overflow/underflow that's inherent + * to determinant computation. + * + * \sa absDeterminant(), MatrixBase::determinant() + */ + typename MatrixType::RealScalar logAbsDeterminant() const; + /** \returns the rank of the matrix of which *this is the QR decomposition. * * \note This is computed at the time of the construction of the QR decomposition. This @@ -238,6 +253,14 @@ typename MatrixType::RealScalar FullPivotingHouseholderQR::absDeterm return ei_abs(m_qr.diagonal().prod()); } +template +typename MatrixType::RealScalar FullPivotingHouseholderQR::logAbsDeterminant() const +{ + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); + return m_qr.diagonal().cwise().abs().cwise().log().sum(); +} + template FullPivotingHouseholderQR& FullPivotingHouseholderQR::compute(const MatrixType& matrix) { diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index d784e0d43..a6430e6f1 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -99,6 +99,7 @@ template void qr_invertible() m1 = m3 * m1 * m3; qr.compute(m1); VERIFY_IS_APPROX(absdet, qr.absDeterminant()); + VERIFY_IS_APPROX(ei_log(absdet), qr.logAbsDeterminant()); } template void qr_verify_assert()