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@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
13namespace Eigen {
14
15namespace 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.
20template<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>
29struct 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
38template<typename T, typename U>
39struct 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  */
60template<typename Derived>
61template<typename OtherDerived>
62typename internal::scalar_product_traits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType
63MatrixBase<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  */
86template<typename Derived>
87template<typename OtherDerived>
88typename internal::traits<Derived>::Scalar
89MatrixBase<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  */
112template<typename Derived>
113EIGEN_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  */
124template<typename Derived>
125inline 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  */
137template<typename Derived>
138inline const typename MatrixBase<Derived>::PlainObject
139MatrixBase<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  */
153template<typename Derived>
154inline void MatrixBase<Derived>::normalize()
155{
156  *this /= norm();
157}
158
159//---------- implementation of other norms ----------
160
161namespace internal {
162
163template<typename Derived, int p>
164struct 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
174template<typename Derived>
175struct 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
183template<typename Derived>
184struct 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
192template<typename Derived>
193struct 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  */
209template<typename Derived>
210template<int p>
211inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
212MatrixBase<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  */
225template<typename Derived>
226template<typename OtherDerived>
227bool 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  */
246template<typename Derived>
247bool 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