add logAbsDeterminant()

move log and exp functors from Array to Core
update documentation

Eigen/src/Array/CwiseOperators.h

@@ -43,38 +43,6 @@ Cwise<ExpressionType>::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<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.

Eigen/src/Array/Functors.h

@@ -69,40 +69,6 @@ struct ei_functor_traits<ei_scalar_sqrt_op<Scalar> >
   };
 };
 
-/** \internal
-  *
-  * \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

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

Eigen/src/Core/Functors.h

@@ -298,6 +298,36 @@ template<typename Scalar>
 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
   *

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
@@ -238,6 +253,14 @@ typename MatrixType::RealScalar FullPivotingHouseholderQR<MatrixType>::absDeterm
   return ei_abs(m_qr.diagonal().prod());
 }
 
+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)
 {

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()