1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2006-2008, 2010 Benoit Jacob <jacob.benoit.1 (at) gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_DOT_H 11 #define EIGEN_DOT_H 12 13 namespace Eigen { 14 15 namespace internal { 16 17 // helper function for dot(). The problem is that if we put that in the body of dot(), then upon calling dot 18 // with mismatched types, the compiler emits errors about failing to instantiate cwiseProduct BEFORE 19 // looking at the static assertions. Thus this is a trick to get better compile errors. 20 template<typename T, typename U, 21 // the NeedToTranspose condition here is taken straight from Assign.h 22 bool NeedToTranspose = T::IsVectorAtCompileTime 23 && U::IsVectorAtCompileTime 24 && ((int(T::RowsAtCompileTime) == 1 && int(U::ColsAtCompileTime) == 1) 25 | // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&". 26 // revert to || as soon as not needed anymore. 27 (int(T::ColsAtCompileTime) == 1 && int(U::RowsAtCompileTime) == 1)) 28 > 29 struct dot_nocheck 30 { 31 typedef typename scalar_product_traits<typename traits<T>::Scalar,typename traits<U>::Scalar>::ReturnType ResScalar; 32 static inline ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b) 33 { 34 return a.template binaryExpr<scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> >(b).sum(); 35 } 36 }; 37 38 template<typename T, typename U> 39 struct dot_nocheck<T, U, true> 40 { 41 typedef typename scalar_product_traits<typename traits<T>::Scalar,typename traits<U>::Scalar>::ReturnType ResScalar; 42 static inline ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b) 43 { 44 return a.transpose().template binaryExpr<scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> >(b).sum(); 45 } 46 }; 47 48 } // end namespace internal 49 50 /** \returns the dot product of *this with other. 51 * 52 * \only_for_vectors 53 * 54 * \note If the scalar type is complex numbers, then this function returns the hermitian 55 * (sesquilinear) dot product, conjugate-linear in the first variable and linear in the 56 * second variable. 57 * 58 * \sa squaredNorm(), norm() 59 */ 60 template<typename Derived> 61 template<typename OtherDerived> 62 typename internal::scalar_product_traits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType 63 MatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const 64 { 65 EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) 66 EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) 67 EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived) 68 typedef internal::scalar_conj_product_op<Scalar,typename OtherDerived::Scalar> func; 69 EIGEN_CHECK_BINARY_COMPATIBILIY(func,Scalar,typename OtherDerived::Scalar); 70 71 eigen_assert(size() == other.size()); 72 73 return internal::dot_nocheck<Derived,OtherDerived>::run(*this, other); 74 } 75 76 #ifdef EIGEN2_SUPPORT 77 /** \returns the dot product of *this with other, with the Eigen2 convention that the dot product is linear in the first variable 78 * (conjugating the second variable). Of course this only makes a difference in the complex case. 79 * 80 * This method is only available in EIGEN2_SUPPORT mode. 81 * 82 * \only_for_vectors 83 * 84 * \sa dot() 85 */ 86 template<typename Derived> 87 template<typename OtherDerived> 88 typename internal::traits<Derived>::Scalar 89 MatrixBase<Derived>::eigen2_dot(const MatrixBase<OtherDerived>& other) const 90 { 91 EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) 92 EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) 93 EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived) 94 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value), 95 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) 96 97 eigen_assert(size() == other.size()); 98 99 return internal::dot_nocheck<OtherDerived,Derived>::run(other,*this); 100 } 101 #endif 102 103 104 //---------- implementation of L2 norm and related functions ---------- 105 106 /** \returns, for vectors, the squared \em l2 norm of \c *this, and for matrices the Frobenius norm. 107 * In both cases, it consists in the sum of the square of all the matrix entries. 108 * For vectors, this is also equals to the dot product of \c *this with itself. 109 * 110 * \sa dot(), norm() 111 */ 112 template<typename Derived> 113 EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::squaredNorm() const 114 { 115 return numext::real((*this).cwiseAbs2().sum()); 116 } 117 118 /** \returns, for vectors, the \em l2 norm of \c *this, and for matrices the Frobenius norm. 119 * In both cases, it consists in the square root of the sum of the square of all the matrix entries. 120 * For vectors, this is also equals to the square root of the dot product of \c *this with itself. 121 * 122 * \sa dot(), squaredNorm() 123 */ 124 template<typename Derived> 125 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::norm() const 126 { 127 using std::sqrt; 128 return sqrt(squaredNorm()); 129 } 130 131 /** \returns an expression of the quotient of *this by its own norm. 132 * 133 * \only_for_vectors 134 * 135 * \sa norm(), normalize() 136 */ 137 template<typename Derived> 138 inline const typename MatrixBase<Derived>::PlainObject 139 MatrixBase<Derived>::normalized() const 140 { 141 typedef typename internal::nested<Derived>::type Nested; 142 typedef typename internal::remove_reference<Nested>::type _Nested; 143 _Nested n(derived()); 144 return n / n.norm(); 145 } 146 147 /** Normalizes the vector, i.e. divides it by its own norm. 148 * 149 * \only_for_vectors 150 * 151 * \sa norm(), normalized() 152 */ 153 template<typename Derived> 154 inline void MatrixBase<Derived>::normalize() 155 { 156 *this /= norm(); 157 } 158 159 //---------- implementation of other norms ---------- 160 161 namespace internal { 162 163 template<typename Derived, int p> 164 struct lpNorm_selector 165 { 166 typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar; 167 static inline RealScalar run(const MatrixBase<Derived>& m) 168 { 169 using std::pow; 170 return pow(m.cwiseAbs().array().pow(p).sum(), RealScalar(1)/p); 171 } 172 }; 173 174 template<typename Derived> 175 struct lpNorm_selector<Derived, 1> 176 { 177 static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m) 178 { 179 return m.cwiseAbs().sum(); 180 } 181 }; 182 183 template<typename Derived> 184 struct lpNorm_selector<Derived, 2> 185 { 186 static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m) 187 { 188 return m.norm(); 189 } 190 }; 191 192 template<typename Derived> 193 struct lpNorm_selector<Derived, Infinity> 194 { 195 static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m) 196 { 197 return m.cwiseAbs().maxCoeff(); 198 } 199 }; 200 201 } // end namespace internal 202 203 /** \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 204 * of the coefficients of *this. If \a p is the special value \a Eigen::Infinity, this function returns the \f$ \ell^\infty \f$ 205 * norm, that is the maximum of the absolute values of the coefficients of *this. 206 * 207 * \sa norm() 208 */ 209 template<typename Derived> 210 template<int p> 211 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real 212 MatrixBase<Derived>::lpNorm() const 213 { 214 return internal::lpNorm_selector<Derived, p>::run(*this); 215 } 216 217 //---------- implementation of isOrthogonal / isUnitary ---------- 218 219 /** \returns true if *this is approximately orthogonal to \a other, 220 * within the precision given by \a prec. 221 * 222 * Example: \include MatrixBase_isOrthogonal.cpp 223 * Output: \verbinclude MatrixBase_isOrthogonal.out 224 */ 225 template<typename Derived> 226 template<typename OtherDerived> 227 bool MatrixBase<Derived>::isOrthogonal 228 (const MatrixBase<OtherDerived>& other, const RealScalar& prec) const 229 { 230 typename internal::nested<Derived,2>::type nested(derived()); 231 typename internal::nested<OtherDerived,2>::type otherNested(other.derived()); 232 return numext::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm(); 233 } 234 235 /** \returns true if *this is approximately an unitary matrix, 236 * within the precision given by \a prec. In the case where the \a Scalar 237 * type is real numbers, a unitary matrix is an orthogonal matrix, whence the name. 238 * 239 * \note This can be used to check whether a family of vectors forms an orthonormal basis. 240 * Indeed, \c m.isUnitary() returns true if and only if the columns (equivalently, the rows) of m form an 241 * orthonormal basis. 242 * 243 * Example: \include MatrixBase_isUnitary.cpp 244 * Output: \verbinclude MatrixBase_isUnitary.out 245 */ 246 template<typename Derived> 247 bool MatrixBase<Derived>::isUnitary(const RealScalar& prec) const 248 { 249 typename Derived::Nested nested(derived()); 250 for(Index i = 0; i < cols(); ++i) 251 { 252 if(!internal::isApprox(nested.col(i).squaredNorm(), static_cast<RealScalar>(1), prec)) 253 return false; 254 for(Index j = 0; j < i; ++j) 255 if(!internal::isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec)) 256 return false; 257 } 258 return true; 259 } 260 261 } // end namespace Eigen 262 263 #endif // EIGEN_DOT_H 264
