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