src/Core/Dot.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2006-2008, 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_DOT_H
0011 #define EIGEN_DOT_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 // helper function for dot(). The problem is that if we put that in the body of dot(), then upon calling dot
0018 // with mismatched types, the compiler emits errors about failing to instantiate cwiseProduct BEFORE
0019 // looking at the static assertions. Thus this is a trick to get better compile errors.
0020 template<typename T, typename U,
0021 // the NeedToTranspose condition here is taken straight from Assign.h
0022          bool NeedToTranspose = T::IsVectorAtCompileTime
0023                 && U::IsVectorAtCompileTime
0024                 && ((int(T::RowsAtCompileTime) == 1 && int(U::ColsAtCompileTime) == 1)
0025                       |  // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&".
0026                          // revert to || as soon as not needed anymore.
0027                     (int(T::ColsAtCompileTime) == 1 && int(U::RowsAtCompileTime) == 1))
0028 >
0029 struct dot_nocheck
0030 {
0031   typedef scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> conj_prod;
0032   typedef typename conj_prod::result_type ResScalar;
0033   EIGEN_DEVICE_FUNC
0034   EIGEN_STRONG_INLINE
0035   static ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b)
0036   {
0037     return a.template binaryExpr<conj_prod>(b).sum();
0038   }
0039 };
0040
0041 template<typename T, typename U>
0042 struct dot_nocheck<T, U, true>
0043 {
0044   typedef scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> conj_prod;
0045   typedef typename conj_prod::result_type ResScalar;
0046   EIGEN_DEVICE_FUNC
0047   EIGEN_STRONG_INLINE
0048   static ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b)
0049   {
0050     return a.transpose().template binaryExpr<conj_prod>(b).sum();
0051   }
0052 };
0053
0054 } // end namespace internal
0055
0056 /** \fn MatrixBase::dot
0057   * \returns the dot product of *this with other.
0058   *
0059   * \only_for_vectors
0060   *
0061   * \note If the scalar type is complex numbers, then this function returns the hermitian
0062   * (sesquilinear) dot product, conjugate-linear in the first variable and linear in the
0063   * second variable.
0064   *
0065   * \sa squaredNorm(), norm()
0066   */
0067 template<typename Derived>
0068 template<typename OtherDerived>
0069 EIGEN_DEVICE_FUNC
0070 EIGEN_STRONG_INLINE
0071 typename ScalarBinaryOpTraits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType
0072 MatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const
0073 {
0074   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
0075   EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
0076   EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
0077 #if !(defined(EIGEN_NO_STATIC_ASSERT) && defined(EIGEN_NO_DEBUG))
0078   typedef internal::scalar_conj_product_op<Scalar,typename OtherDerived::Scalar> func;
0079   EIGEN_CHECK_BINARY_COMPATIBILIY(func,Scalar,typename OtherDerived::Scalar);
0080 #endif
0081
0082   eigen_assert(size() == other.size());
0083
0084   return internal::dot_nocheck<Derived,OtherDerived>::run(*this, other);
0085 }
0086
0087 //---------- implementation of L2 norm and related functions ----------
0088
0089 /** \returns, for vectors, the squared \em l2 norm of \c *this, and for matrices the squared Frobenius norm.
0090   * In both cases, it consists in the sum of the square of all the matrix entries.
0091   * For vectors, this is also equals to the dot product of \c *this with itself.
0092   *
0093   * \sa dot(), norm(), lpNorm()
0094   */
0095 template<typename Derived>
0096 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::squaredNorm() const
0097 {
0098   return numext::real((*this).cwiseAbs2().sum());
0099 }
0100
0101 /** \returns, for vectors, the \em l2 norm of \c *this, and for matrices the Frobenius norm.
0102   * In both cases, it consists in the square root of the sum of the square of all the matrix entries.
0103   * For vectors, this is also equals to the square root of the dot product of \c *this with itself.
0104   *
0105   * \sa lpNorm(), dot(), squaredNorm()
0106   */
0107 template<typename Derived>
0108 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::norm() const
0109 {
0110   return numext::sqrt(squaredNorm());
0111 }
0112
0113 /** \returns an expression of the quotient of \c *this by its own norm.
0114   *
0115   * \warning If the input vector is too small (i.e., this->norm()==0),
0116   *          then this function returns a copy of the input.
0117   *
0118   * \only_for_vectors
0119   *
0120   * \sa norm(), normalize()
0121   */
0122 template<typename Derived>
0123 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::PlainObject
0124 MatrixBase<Derived>::normalized() const
0125 {
0126   typedef typename internal::nested_eval<Derived,2>::type _Nested;
0127   _Nested n(derived());
0128   RealScalar z = n.squaredNorm();
0129   // NOTE: after extensive benchmarking, this conditional does not impact performance, at least on recent x86 CPU
0130   if(z>RealScalar(0))
0131     return n / numext::sqrt(z);
0132   else
0133     return n;
0134 }
0135
0136 /** Normalizes the vector, i.e. divides it by its own norm.
0137   *
0138   * \only_for_vectors
0139   *
0140   * \warning If the input vector is too small (i.e., this->norm()==0), then \c *this is left unchanged.
0141   *
0142   * \sa norm(), normalized()
0143   */
0144 template<typename Derived>
0145 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase<Derived>::normalize()
0146 {
0147   RealScalar z = squaredNorm();
0148   // NOTE: after extensive benchmarking, this conditional does not impact performance, at least on recent x86 CPU
0149   if(z>RealScalar(0))
0150     derived() /= numext::sqrt(z);
0151 }
0152
0153 /** \returns an expression of the quotient of \c *this by its own norm while avoiding underflow and overflow.
0154   *
0155   * \only_for_vectors
0156   *
0157   * This method is analogue to the normalized() method, but it reduces the risk of
0158   * underflow and overflow when computing the norm.
0159   *
0160   * \warning If the input vector is too small (i.e., this->norm()==0),
0161   *          then this function returns a copy of the input.
0162   *
0163   * \sa stableNorm(), stableNormalize(), normalized()
0164   */
0165 template<typename Derived>
0166 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::PlainObject
0167 MatrixBase<Derived>::stableNormalized() const
0168 {
0169   typedef typename internal::nested_eval<Derived,3>::type _Nested;
0170   _Nested n(derived());
0171   RealScalar w = n.cwiseAbs().maxCoeff();
0172   RealScalar z = (n/w).squaredNorm();
0173   if(z>RealScalar(0))
0174     return n / (numext::sqrt(z)*w);
0175   else
0176     return n;
0177 }
0178
0179 /** Normalizes the vector while avoid underflow and overflow
0180   *
0181   * \only_for_vectors
0182   *
0183   * This method is analogue to the normalize() method, but it reduces the risk of
0184   * underflow and overflow when computing the norm.
0185   *
0186   * \warning If the input vector is too small (i.e., this->norm()==0), then \c *this is left unchanged.
0187   *
0188   * \sa stableNorm(), stableNormalized(), normalize()
0189   */
0190 template<typename Derived>
0191 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase<Derived>::stableNormalize()
0192 {
0193   RealScalar w = cwiseAbs().maxCoeff();
0194   RealScalar z = (derived()/w).squaredNorm();
0195   if(z>RealScalar(0))
0196     derived() /= numext::sqrt(z)*w;
0197 }
0198
0199 //---------- implementation of other norms ----------
0200
0201 namespace internal {
0202
0203 template<typename Derived, int p>
0204 struct lpNorm_selector
0205 {
0206   typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
0207   EIGEN_DEVICE_FUNC
0208   static inline RealScalar run(const MatrixBase<Derived>& m)
0209   {
0210     EIGEN_USING_STD(pow)
0211     return pow(m.cwiseAbs().array().pow(p).sum(), RealScalar(1)/p);
0212   }
0213 };
0214
0215 template<typename Derived>
0216 struct lpNorm_selector<Derived, 1>
0217 {
0218   EIGEN_DEVICE_FUNC
0219   static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
0220   {
0221     return m.cwiseAbs().sum();
0222   }
0223 };
0224
0225 template<typename Derived>
0226 struct lpNorm_selector<Derived, 2>
0227 {
0228   EIGEN_DEVICE_FUNC
0229   static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
0230   {
0231     return m.norm();
0232   }
0233 };
0234
0235 template<typename Derived>
0236 struct lpNorm_selector<Derived, Infinity>
0237 {
0238   typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
0239   EIGEN_DEVICE_FUNC
0240   static inline RealScalar run(const MatrixBase<Derived>& m)
0241   {
0242     if(Derived::SizeAtCompileTime==0 || (Derived::SizeAtCompileTime==Dynamic && m.size()==0))
0243       return RealScalar(0);
0244     return m.cwiseAbs().maxCoeff();
0245   }
0246 };
0247
0248 } // end namespace internal
0249
0250 /** \returns the \b coefficient-wise \f$ \ell^p \f$ norm of \c *this, that is, returns the p-th root of the sum of the p-th powers of the absolute values
0251   *          of the coefficients of \c *this. If \a p is the special value \a Eigen::Infinity, this function returns the \f$ \ell^\infty \f$
0252   *          norm, that is the maximum of the absolute values of the coefficients of \c *this.
0253   *
0254   * In all cases, if \c *this is empty, then the value 0 is returned.
0255   *
0256   * \note For matrices, this function does not compute the <a href="https://en.wikipedia.org/wiki/Operator_norm">operator-norm</a>. That is, if \c *this is a matrix, then its coefficients are interpreted as a 1D vector. Nonetheless, you can easily compute the 1-norm and \f$\infty\f$-norm matrix operator norms using \link TutorialReductionsVisitorsBroadcastingReductionsNorm partial reductions \endlink.
0257   *
0258   * \sa norm()
0259   */
0260 template<typename Derived>
0261 template<int p>
0262 #ifndef EIGEN_PARSED_BY_DOXYGEN
0263 EIGEN_DEVICE_FUNC inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
0264 #else
0265 EIGEN_DEVICE_FUNC MatrixBase<Derived>::RealScalar
0266 #endif
0267 MatrixBase<Derived>::lpNorm() const
0268 {
0269   return internal::lpNorm_selector<Derived, p>::run(*this);
0270 }
0271
0272 //---------- implementation of isOrthogonal / isUnitary ----------
0273
0274 /** \returns true if *this is approximately orthogonal to \a other,
0275   *          within the precision given by \a prec.
0276   *
0277   * Example: \include MatrixBase_isOrthogonal.cpp
0278   * Output: \verbinclude MatrixBase_isOrthogonal.out
0279   */
0280 template<typename Derived>
0281 template<typename OtherDerived>
0282 bool MatrixBase<Derived>::isOrthogonal
0283 (const MatrixBase<OtherDerived>& other, const RealScalar& prec) const
0284 {
0285   typename internal::nested_eval<Derived,2>::type nested(derived());
0286   typename internal::nested_eval<OtherDerived,2>::type otherNested(other.derived());
0287   return numext::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm();
0288 }
0289
0290 /** \returns true if *this is approximately an unitary matrix,
0291   *          within the precision given by \a prec. In the case where the \a Scalar
0292   *          type is real numbers, a unitary matrix is an orthogonal matrix, whence the name.
0293   *
0294   * \note This can be used to check whether a family of vectors forms an orthonormal basis.
0295   *       Indeed, \c m.isUnitary() returns true if and only if the columns (equivalently, the rows) of m form an
0296   *       orthonormal basis.
0297   *
0298   * Example: \include MatrixBase_isUnitary.cpp
0299   * Output: \verbinclude MatrixBase_isUnitary.out
0300   */
0301 template<typename Derived>
0302 bool MatrixBase<Derived>::isUnitary(const RealScalar& prec) const
0303 {
0304   typename internal::nested_eval<Derived,1>::type self(derived());
0305   for(Index i = 0; i < cols(); ++i)
0306   {
0307     if(!internal::isApprox(self.col(i).squaredNorm(), static_cast<RealScalar>(1), prec))
0308       return false;
0309     for(Index j = 0; j < i; ++j)
0310       if(!internal::isMuchSmallerThan(self.col(i).dot(self.col(j)), static_cast<Scalar>(1), prec))
0311         return false;
0312   }
0313   return true;
0314 }
0315
0316 } // end namespace Eigen
0317
0318 #endif // EIGEN_DOT_H