Dot.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2006-2008, 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_DOT_H
00011 #define EIGEN_DOT_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 // helper function for dot(). The problem is that if we put that in the body of dot(), then upon calling dot
00018 // with mismatched types, the compiler emits errors about failing to instantiate cwiseProduct BEFORE
00019 // looking at the static assertions. Thus this is a trick to get better compile errors.
00020 template<typename T, typename U,
00021 // the NeedToTranspose condition here is taken straight from Assign.h
00022          bool NeedToTranspose = T::IsVectorAtCompileTime
00023                 && U::IsVectorAtCompileTime
00024                 && ((int(T::RowsAtCompileTime) == 1 && int(U::ColsAtCompileTime) == 1)
00025                       |  // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&".
00026                          // revert to || as soon as not needed anymore.
00027                     (int(T::ColsAtCompileTime) == 1 && int(U::RowsAtCompileTime) == 1))
00028 >
00029 struct dot_nocheck
00030 {
00031   typedef typename scalar_product_traits<typename traits<T>::Scalar,typename traits<U>::Scalar>::ReturnType ResScalar;
00032   static inline ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b)
00033   {
00034     return a.template binaryExpr<scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> >(b).sum();
00035   }
00036 };
00037 
00038 template<typename T, typename U>
00039 struct dot_nocheck<T, U, true>
00040 {
00041   typedef typename scalar_product_traits<typename traits<T>::Scalar,typename traits<U>::Scalar>::ReturnType ResScalar;
00042   static inline ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b)
00043   {
00044     return a.transpose().template binaryExpr<scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> >(b).sum();
00045   }
00046 };
00047 
00048 } // end namespace internal
00049 
00050 /** \returns the dot product of *this with other.
00051   *
00052   * \only_for_vectors
00053   *
00054   * \note If the scalar type is complex numbers, then this function returns the hermitian
00055   * (sesquilinear) dot product, conjugate-linear in the first variable and linear in the
00056   * second variable.
00057   *
00058   * \sa squaredNorm(), norm()
00059   */
00060 template<typename Derived>
00061 template<typename OtherDerived>
00062 typename internal::scalar_product_traits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType
00063 MatrixBase<Derived>::dot (const MatrixBase<OtherDerived>& other) const
00064 {
00065   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
00066   EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
00067   EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
00068   typedef internal::scalar_conj_product_op<Scalar,typename OtherDerived::Scalar> func;
00069   EIGEN_CHECK_BINARY_COMPATIBILIY(func,Scalar,typename OtherDerived::Scalar);
00070 
00071   eigen_assert(size() == other.size());
00072 
00073   return internal::dot_nocheck<Derived,OtherDerived>::run(*this, other);
00074 }
00075 
00076 #ifdef EIGEN2_SUPPORT
00077 /** \returns the dot product of *this with other, with the Eigen2 convention that the dot product is linear in the first variable
00078   * (conjugating the second variable). Of course this only makes a difference in the complex case.
00079   *
00080   * This method is only available in EIGEN2_SUPPORT mode.
00081   *
00082   * \only_for_vectors
00083   *
00084   * \sa dot()
00085   */
00086 template<typename Derived>
00087 template<typename OtherDerived>
00088 typename internal::traits<Derived>::Scalar
00089 MatrixBase<Derived>::eigen2_dot (const MatrixBase<OtherDerived>& other) const
00090 {
00091   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
00092   EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
00093   EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
00094   EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value),
00095     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00096 
00097   eigen_assert(size() == other.size());
00098 
00099   return internal::dot_nocheck<OtherDerived,Derived>::run(other,*this);
00100 }
00101 #endif
00102 
00103 
00104 //---------- implementation of L2 norm and related functions ----------
00105 
00106 /** \returns, for vectors, the squared \em l2 norm of \c *this, and for matrices the Frobenius norm.
00107   * In both cases, it consists in the sum of the square of all the matrix entries.
00108   * For vectors, this is also equals to the dot product of \c *this with itself.
00109   *
00110   * \sa dot(), norm()
00111   */
00112 template<typename Derived>
00113 EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::squaredNorm () const
00114 {
00115   return numext::real((*this).cwiseAbs2().sum());
00116 }
00117 
00118 /** \returns, for vectors, the \em l2 norm of \c *this, and for matrices the Frobenius norm.
00119   * In both cases, it consists in the square root of the sum of the square of all the matrix entries.
00120   * For vectors, this is also equals to the square root of the dot product of \c *this with itself.
00121   *
00122   * \sa dot(), squaredNorm()
00123   */
00124 template<typename Derived>
00125 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::norm () const
00126 {
00127   using std::sqrt;
00128   return sqrt(squaredNorm());
00129 }
00130 
00131 /** \returns an expression of the quotient of *this by its own norm.
00132   *
00133   * \only_for_vectors
00134   *
00135   * \sa norm(), normalize()
00136   */
00137 template<typename Derived>
00138 inline const typename MatrixBase<Derived>::PlainObject
00139 MatrixBase<Derived>::normalized () const
00140 {
00141   typedef typename internal::nested<Derived>::type Nested;
00142   typedef typename internal::remove_reference<Nested>::type _Nested;
00143   _Nested n(derived());
00144   return n / n.norm();
00145 }
00146 
00147 /** Normalizes the vector, i.e. divides it by its own norm.
00148   *
00149   * \only_for_vectors
00150   *
00151   * \sa norm(), normalized()
00152   */
00153 template<typename Derived>
00154 inline void MatrixBase<Derived>::normalize()
00155 {
00156   *this /= norm();
00157 }
00158 
00159 //---------- implementation of other norms ----------
00160 
00161 namespace internal {
00162 
00163 template<typename Derived, int p>
00164 struct lpNorm_selector
00165 {
00166   typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
00167   static inline RealScalar run(const MatrixBase<Derived>& m)
00168   {
00169     using std::pow;
00170     return pow(m.cwiseAbs().array ().pow(p).sum(), RealScalar(1)/p);
00171   }
00172 };
00173 
00174 template<typename Derived>
00175 struct lpNorm_selector<Derived, 1>
00176 {
00177   static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
00178   {
00179     return m.cwiseAbs().sum();
00180   }
00181 };
00182 
00183 template<typename Derived>
00184 struct lpNorm_selector<Derived, 2>
00185 {
00186   static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
00187   {
00188     return m.norm();
00189   }
00190 };
00191 
00192 template<typename Derived>
00193 struct lpNorm_selector<Derived, Infinity>
00194 {
00195   static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
00196   {
00197     return m.cwiseAbs().maxCoeff();
00198   }
00199 };
00200 
00201 } // end namespace internal
00202 
00203 /** \returns the \f$ \ell^p \f$ norm of *this, that is, returns the p-th root of the sum of the p-th powers of the absolute values
00204   *          of the coefficients of *this. If \a p is the special value \a Eigen::Infinity, this function returns the \f$ \ell^\infty \f$
00205   *          norm, that is the maximum of the absolute values of the coefficients of *this.
00206   *
00207   * \sa norm()
00208   */
00209 template<typename Derived>
00210 template<int p>
00211 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
00212 MatrixBase<Derived>::lpNorm () const
00213 {
00214   return internal::lpNorm_selector<Derived, p>::run(*this);
00215 }
00216 
00217 //---------- implementation of isOrthogonal / isUnitary ----------
00218 
00219 /** \returns true if *this is approximately orthogonal to \a other,
00220   *          within the precision given by \a prec.
00221   *
00222   * Example: \include MatrixBase_isOrthogonal.cpp
00223   * Output: \verbinclude MatrixBase_isOrthogonal.out
00224   */
00225 template<typename Derived>
00226 template<typename OtherDerived>
00227 bool MatrixBase<Derived>::isOrthogonal 
00228 (const MatrixBase<OtherDerived>& other, const RealScalar& prec) const
00229 {
00230   typename internal::nested<Derived,2>::type nested(derived());
00231   typename internal::nested<OtherDerived,2>::type otherNested(other.derived());
00232   return numext::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm();
00233 }
00234 
00235 /** \returns true if *this is approximately an unitary matrix,
00236   *          within the precision given by \a prec. In the case where the \a Scalar
00237   *          type is real numbers, a unitary matrix is an orthogonal matrix, whence the name.
00238   *
00239   * \note This can be used to check whether a family of vectors forms an orthonormal basis.
00240   *       Indeed, \c m.isUnitary() returns true if and only if the columns (equivalently, the rows) of m form an
00241   *       orthonormal basis.
00242   *
00243   * Example: \include MatrixBase_isUnitary.cpp
00244   * Output: \verbinclude MatrixBase_isUnitary.out
00245   */
00246 template<typename Derived>
00247 bool MatrixBase<Derived>::isUnitary (const RealScalar& prec) const
00248 {
00249   typename Derived::Nested nested(derived());
00250   for(Index i = 0; i < cols(); ++i)
00251   {
00252     if(!internal::isApprox(nested.col(i).squaredNorm(), static_cast<RealScalar>(1), prec))
00253       return false;
00254     for(Index j = 0; j < i; ++j)
00255       if(!internal::isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec))
00256         return false;
00257   }
00258   return true;
00259 }
00260 
00261 } // end namespace Eigen
00262 
00263 #endif // EIGEN_DOT_H