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 scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> conj_prod;
32  typedef typename conj_prod::result_type ResScalar;
33  EIGEN_DEVICE_FUNC
34  static inline ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b)
35  {
36    return a.template binaryExpr<conj_prod>(b).sum();
37  }
38};
39
40template<typename T, typename U>
41struct dot_nocheck<T, U, true>
42{
43  typedef scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> conj_prod;
44  typedef typename conj_prod::result_type ResScalar;
45  EIGEN_DEVICE_FUNC
46  static inline ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b)
47  {
48    return a.transpose().template binaryExpr<conj_prod>(b).sum();
49  }
50};
51
52} // end namespace internal
53
54/** \fn MatrixBase::dot
55  * \returns the dot product of *this with other.
56  *
57  * \only_for_vectors
58  *
59  * \note If the scalar type is complex numbers, then this function returns the hermitian
60  * (sesquilinear) dot product, conjugate-linear in the first variable and linear in the
61  * second variable.
62  *
63  * \sa squaredNorm(), norm()
64  */
65template<typename Derived>
66template<typename OtherDerived>
67EIGEN_DEVICE_FUNC
68typename ScalarBinaryOpTraits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType
69MatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const
70{
71  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
72  EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
73  EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
74#if !(defined(EIGEN_NO_STATIC_ASSERT) && defined(EIGEN_NO_DEBUG))
75  typedef internal::scalar_conj_product_op<Scalar,typename OtherDerived::Scalar> func;
76  EIGEN_CHECK_BINARY_COMPATIBILIY(func,Scalar,typename OtherDerived::Scalar);
77#endif
78
79  eigen_assert(size() == other.size());
80
81  return internal::dot_nocheck<Derived,OtherDerived>::run(*this, other);
82}
83
84//---------- implementation of L2 norm and related functions ----------
85
86/** \returns, for vectors, the squared \em l2 norm of \c *this, and for matrices the Frobenius norm.
87  * In both cases, it consists in the sum of the square of all the matrix entries.
88  * For vectors, this is also equals to the dot product of \c *this with itself.
89  *
90  * \sa dot(), norm(), lpNorm()
91  */
92template<typename Derived>
93EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::squaredNorm() const
94{
95  return numext::real((*this).cwiseAbs2().sum());
96}
97
98/** \returns, for vectors, the \em l2 norm of \c *this, and for matrices the Frobenius norm.
99  * In both cases, it consists in the square root of the sum of the square of all the matrix entries.
100  * For vectors, this is also equals to the square root of the dot product of \c *this with itself.
101  *
102  * \sa lpNorm(), dot(), squaredNorm()
103  */
104template<typename Derived>
105inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::norm() const
106{
107  return numext::sqrt(squaredNorm());
108}
109
110/** \returns an expression of the quotient of \c *this by its own norm.
111  *
112  * \warning If the input vector is too small (i.e., this->norm()==0),
113  *          then this function returns a copy of the input.
114  *
115  * \only_for_vectors
116  *
117  * \sa norm(), normalize()
118  */
119template<typename Derived>
120inline const typename MatrixBase<Derived>::PlainObject
121MatrixBase<Derived>::normalized() const
122{
123  typedef typename internal::nested_eval<Derived,2>::type _Nested;
124  _Nested n(derived());
125  RealScalar z = n.squaredNorm();
126  // NOTE: after extensive benchmarking, this conditional does not impact performance, at least on recent x86 CPU
127  if(z>RealScalar(0))
128    return n / numext::sqrt(z);
129  else
130    return n;
131}
132
133/** Normalizes the vector, i.e. divides it by its own norm.
134  *
135  * \only_for_vectors
136  *
137  * \warning If the input vector is too small (i.e., this->norm()==0), then \c *this is left unchanged.
138  *
139  * \sa norm(), normalized()
140  */
141template<typename Derived>
142inline void MatrixBase<Derived>::normalize()
143{
144  RealScalar z = squaredNorm();
145  // NOTE: after extensive benchmarking, this conditional does not impact performance, at least on recent x86 CPU
146  if(z>RealScalar(0))
147    derived() /= numext::sqrt(z);
148}
149
150/** \returns an expression of the quotient of \c *this by its own norm while avoiding underflow and overflow.
151  *
152  * \only_for_vectors
153  *
154  * This method is analogue to the normalized() method, but it reduces the risk of
155  * underflow and overflow when computing the norm.
156  *
157  * \warning If the input vector is too small (i.e., this->norm()==0),
158  *          then this function returns a copy of the input.
159  *
160  * \sa stableNorm(), stableNormalize(), normalized()
161  */
162template<typename Derived>
163inline const typename MatrixBase<Derived>::PlainObject
164MatrixBase<Derived>::stableNormalized() const
165{
166  typedef typename internal::nested_eval<Derived,3>::type _Nested;
167  _Nested n(derived());
168  RealScalar w = n.cwiseAbs().maxCoeff();
169  RealScalar z = (n/w).squaredNorm();
170  if(z>RealScalar(0))
171    return n / (numext::sqrt(z)*w);
172  else
173    return n;
174}
175
176/** Normalizes the vector while avoid underflow and overflow
177  *
178  * \only_for_vectors
179  *
180  * This method is analogue to the normalize() method, but it reduces the risk of
181  * underflow and overflow when computing the norm.
182  *
183  * \warning If the input vector is too small (i.e., this->norm()==0), then \c *this is left unchanged.
184  *
185  * \sa stableNorm(), stableNormalized(), normalize()
186  */
187template<typename Derived>
188inline void MatrixBase<Derived>::stableNormalize()
189{
190  RealScalar w = cwiseAbs().maxCoeff();
191  RealScalar z = (derived()/w).squaredNorm();
192  if(z>RealScalar(0))
193    derived() /= numext::sqrt(z)*w;
194}
195
196//---------- implementation of other norms ----------
197
198namespace internal {
199
200template<typename Derived, int p>
201struct lpNorm_selector
202{
203  typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
204  EIGEN_DEVICE_FUNC
205  static inline RealScalar run(const MatrixBase<Derived>& m)
206  {
207    EIGEN_USING_STD_MATH(pow)
208    return pow(m.cwiseAbs().array().pow(p).sum(), RealScalar(1)/p);
209  }
210};
211
212template<typename Derived>
213struct lpNorm_selector<Derived, 1>
214{
215  EIGEN_DEVICE_FUNC
216  static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
217  {
218    return m.cwiseAbs().sum();
219  }
220};
221
222template<typename Derived>
223struct lpNorm_selector<Derived, 2>
224{
225  EIGEN_DEVICE_FUNC
226  static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
227  {
228    return m.norm();
229  }
230};
231
232template<typename Derived>
233struct lpNorm_selector<Derived, Infinity>
234{
235  typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
236  EIGEN_DEVICE_FUNC
237  static inline RealScalar run(const MatrixBase<Derived>& m)
238  {
239    if(Derived::SizeAtCompileTime==0 || (Derived::SizeAtCompileTime==Dynamic && m.size()==0))
240      return RealScalar(0);
241    return m.cwiseAbs().maxCoeff();
242  }
243};
244
245} // end namespace internal
246
247/** \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
248  *          of the coefficients of \c *this. If \a p is the special value \a Eigen::Infinity, this function returns the \f$ \ell^\infty \f$
249  *          norm, that is the maximum of the absolute values of the coefficients of \c *this.
250  *
251  * In all cases, if \c *this is empty, then the value 0 is returned.
252  *
253  * \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.
254  *
255  * \sa norm()
256  */
257template<typename Derived>
258template<int p>
259#ifndef EIGEN_PARSED_BY_DOXYGEN
260inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
261#else
262MatrixBase<Derived>::RealScalar
263#endif
264MatrixBase<Derived>::lpNorm() const
265{
266  return internal::lpNorm_selector<Derived, p>::run(*this);
267}
268
269//---------- implementation of isOrthogonal / isUnitary ----------
270
271/** \returns true if *this is approximately orthogonal to \a other,
272  *          within the precision given by \a prec.
273  *
274  * Example: \include MatrixBase_isOrthogonal.cpp
275  * Output: \verbinclude MatrixBase_isOrthogonal.out
276  */
277template<typename Derived>
278template<typename OtherDerived>
279bool MatrixBase<Derived>::isOrthogonal
280(const MatrixBase<OtherDerived>& other, const RealScalar& prec) const
281{
282  typename internal::nested_eval<Derived,2>::type nested(derived());
283  typename internal::nested_eval<OtherDerived,2>::type otherNested(other.derived());
284  return numext::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm();
285}
286
287/** \returns true if *this is approximately an unitary matrix,
288  *          within the precision given by \a prec. In the case where the \a Scalar
289  *          type is real numbers, a unitary matrix is an orthogonal matrix, whence the name.
290  *
291  * \note This can be used to check whether a family of vectors forms an orthonormal basis.
292  *       Indeed, \c m.isUnitary() returns true if and only if the columns (equivalently, the rows) of m form an
293  *       orthonormal basis.
294  *
295  * Example: \include MatrixBase_isUnitary.cpp
296  * Output: \verbinclude MatrixBase_isUnitary.out
297  */
298template<typename Derived>
299bool MatrixBase<Derived>::isUnitary(const RealScalar& prec) const
300{
301  typename internal::nested_eval<Derived,1>::type self(derived());
302  for(Index i = 0; i < cols(); ++i)
303  {
304    if(!internal::isApprox(self.col(i).squaredNorm(), static_cast<RealScalar>(1), prec))
305      return false;
306    for(Index j = 0; j < i; ++j)
307      if(!internal::isMuchSmallerThan(self.col(i).dot(self.col(j)), static_cast<Scalar>(1), prec))
308        return false;
309  }
310  return true;
311}
312
313} // end namespace Eigen
314
315#endif // EIGEN_DOT_H
316