diff options
| -rw-r--r-- | Eigen/src/Array/CwiseOperators.h | 32 | ||||
| -rw-r--r-- | Eigen/src/Array/Functors.h | 34 | ||||
| -rw-r--r-- | Eigen/src/Core/CwiseUnaryOp.h | 29 | ||||
| -rw-r--r-- | Eigen/src/Core/Functors.h | 30 | ||||
| -rw-r--r-- | Eigen/src/QR/FullPivotingHouseholderQR.h | 31 | ||||
| -rw-r--r-- | test/qr_fullpivoting.cpp | 1 |
6 files changed, 87 insertions, 70 deletions
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 @@ -45,38 +45,6 @@ Cwise<ExpressionType>::sqrt() const /** \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<typename ExpressionType> -inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_exp_op) -Cwise<ExpressionType>::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<typename ExpressionType> -inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_log_op) -Cwise<ExpressionType>::log() const -{ - return _expression(); -} - -/** \array_module - * * \returns an expression of the coefficient-wise cosine of *this. * * Example: \include Cwise_cos.cpp 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 @@ -73,40 +73,6 @@ struct ei_functor_traits<ei_scalar_sqrt_op<Scalar> > * * \array_module * - * \brief Template functor to compute the exponential of a scalar - * - * \sa class CwiseUnaryOp, Cwise::exp() - */ -template<typename Scalar> 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<Scalar>::type Packet; - inline Packet packetOp(const Packet& a) const { return ei_pexp(a); } -}; -template<typename Scalar> -struct ei_functor_traits<ei_scalar_exp_op<Scalar> > -{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::HasExp }; }; - -/** \internal - * - * \array_module - * - * \brief Template functor to compute the logarithm of a scalar - * - * \sa class CwiseUnaryOp, Cwise::log() - */ -template<typename Scalar> 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<Scalar>::type Packet; - inline Packet packetOp(const Packet& a) const { return ei_plog(a); } -}; -template<typename Scalar> -struct ei_functor_traits<ei_scalar_log_op<Scalar> > -{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::HasLog }; }; - -/** \internal - * - * \array_module - * * \brief Template functor to compute the cosine of a scalar * * \sa class CwiseUnaryOp, Cwise::cos() 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<Derived>::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<typename ExpressionType> +inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_exp_op) +Cwise<ExpressionType>::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<typename ExpressionType> +inline const EIGEN_CWISE_UNOP_RETURN_TYPE(ei_scalar_log_op) +Cwise<ExpressionType>::log() const +{ + return _expression(); +} + + /** \relates MatrixBase */ template<typename Derived> EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::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 @@ -299,6 +299,36 @@ struct ei_functor_traits<ei_scalar_imag_op<Scalar> > { enum { Cost = 0, PacketAccess = false }; }; /** \internal + * + * \brief Template functor to compute the exponential of a scalar + * + * \sa class CwiseUnaryOp, Cwise::exp() + */ +template<typename Scalar> 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<Scalar>::type Packet; + inline Packet packetOp(const Packet& a) const { return ei_pexp(a); } +}; +template<typename Scalar> +struct ei_functor_traits<ei_scalar_exp_op<Scalar> > +{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::HasExp }; }; + +/** \internal + * + * \brief Template functor to compute the logarithm of a scalar + * + * \sa class CwiseUnaryOp, Cwise::log() + */ +template<typename Scalar> 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<Scalar>::type Packet; + inline Packet packetOp(const Packet& a) const { return ei_plog(a); } +}; +template<typename Scalar> +struct ei_functor_traits<ei_scalar_log_op<Scalar> > +{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::HasLog }; }; + +/** \internal * \brief Template functor to multiply a scalar by a fixed other one * * \sa class CwiseUnaryOp, MatrixBase::operator*, MatrixBase::operator/ 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<typename MatrixType> class FullPivotingHouseholderQR { @@ -125,11 +125,26 @@ template<typename MatrixType> 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 @@ -239,6 +254,14 @@ typename MatrixType::RealScalar FullPivotingHouseholderQR<MatrixType>::absDeterm } template<typename MatrixType> +typename MatrixType::RealScalar FullPivotingHouseholderQR<MatrixType>::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<typename MatrixType> FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::compute(const MatrixType& matrix) { int rows = matrix.rows(); 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<typename MatrixType> 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<typename MatrixType> void qr_verify_assert() |