Dot.h
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 
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 
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
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 
86 template<typename Derived>
87 template<typename OtherDerived>
88 typename internal::traits<Derived>::Scalar
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 
112 template<typename Derived>
114 {
115  return internal::real((*this).cwiseAbs2().sum());
116 }
117 
124 template<typename Derived>
126 {
127  return internal::sqrt(squaredNorm());
128 }
129 
136 template<typename Derived>
137 inline const typename MatrixBase<Derived>::PlainObject
139 {
140  typedef typename internal::nested<Derived>::type Nested;
141  typedef typename internal::remove_reference<Nested>::type _Nested;
142  _Nested n(derived());
143  return n / n.norm();
144 }
145 
152 template<typename Derived>
154 {
155  *this /= norm();
156 }
157 
158 //---------- implementation of other norms ----------
159 
160 namespace internal {
161 
162 template<typename Derived, int p>
163 struct lpNorm_selector
164 {
165  typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
166  static inline RealScalar run(const MatrixBase<Derived>& m)
167  {
168  return pow(m.cwiseAbs().array().pow(p).sum(), RealScalar(1)/p);
169  }
170 };
171 
172 template<typename Derived>
173 struct lpNorm_selector<Derived, 1>
174 {
175  static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
176  {
177  return m.cwiseAbs().sum();
178  }
179 };
180 
181 template<typename Derived>
182 struct lpNorm_selector<Derived, 2>
183 {
184  static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
185  {
186  return m.norm();
187  }
188 };
189 
190 template<typename Derived>
191 struct lpNorm_selector<Derived, Infinity>
192 {
193  static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
194  {
195  return m.cwiseAbs().maxCoeff();
196  }
197 };
198 
199 } // end namespace internal
200 
207 template<typename Derived>
208 template<int p>
209 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
211 {
212  return internal::lpNorm_selector<Derived, p>::run(*this);
213 }
214 
215 //---------- implementation of isOrthogonal / isUnitary ----------
216 
223 template<typename Derived>
224 template<typename OtherDerived>
226 (const MatrixBase<OtherDerived>& other, RealScalar prec) const
227 {
228  typename internal::nested<Derived,2>::type nested(derived());
229  typename internal::nested<OtherDerived,2>::type otherNested(other.derived());
230  return internal::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm();
231 }
232 
244 template<typename Derived>
245 bool MatrixBase<Derived>::isUnitary(RealScalar prec) const
246 {
247  typename Derived::Nested nested(derived());
248  for(Index i = 0; i < cols(); ++i)
249  {
250  if(!internal::isApprox(nested.col(i).squaredNorm(), static_cast<RealScalar>(1), prec))
251  return false;
252  for(Index j = 0; j < i; ++j)
253  if(!internal::isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec))
254  return false;
255  }
256  return true;
257 }
258 
259 } // end namespace Eigen
260 
261 #endif // EIGEN_DOT_H