Michael Ernst Peter / Eigen

src/Core/Dot.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

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