Dot.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2006-2008, 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_DOT_H
00026 #define EIGEN_DOT_H
00027 
00028 namespace Eigen { 
00029 
00030 namespace internal {
00031 
00032 // helper function for dot(). The problem is that if we put that in the body of dot(), then upon calling dot
00033 // with mismatched types, the compiler emits errors about failing to instantiate cwiseProduct BEFORE
00034 // looking at the static assertions. Thus this is a trick to get better compile errors.
00035 template<typename T, typename U,
00036 // the NeedToTranspose condition here is taken straight from Assign.h
00037          bool NeedToTranspose = T::IsVectorAtCompileTime
00038                 && U::IsVectorAtCompileTime
00039                 && ((int(T::RowsAtCompileTime) == 1 && int(U::ColsAtCompileTime) == 1)
00040                       |  // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&".
00041                          // revert to || as soon as not needed anymore.
00042                     (int(T::ColsAtCompileTime) == 1 && int(U::RowsAtCompileTime) == 1))
00043 >
00044 struct dot_nocheck
00045 {
00046   typedef typename scalar_product_traits<typename traits<T>::Scalar,typename traits<U>::Scalar>::ReturnType ResScalar;
00047   static inline ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b)
00048   {
00049     return a.template binaryExpr<scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> >(b).sum();
00050   }
00051 };
00052 
00053 template<typename T, typename U>
00054 struct dot_nocheck<T, U, true>
00055 {
00056   typedef typename scalar_product_traits<typename traits<T>::Scalar,typename traits<U>::Scalar>::ReturnType ResScalar;
00057   static inline ResScalar run(const MatrixBase<T>& a, const MatrixBase<U>& b)
00058   {
00059     return a.transpose().template binaryExpr<scalar_conj_product_op<typename traits<T>::Scalar,typename traits<U>::Scalar> >(b).sum();
00060   }
00061 };
00062 
00063 } // end namespace internal
00064 
00075 template<typename Derived>
00076 template<typename OtherDerived>
00077 typename internal::scalar_product_traits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType
00078 MatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const
00079 {
00080   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
00081   EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
00082   EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
00083   typedef internal::scalar_conj_product_op<Scalar,typename OtherDerived::Scalar> func;
00084   EIGEN_CHECK_BINARY_COMPATIBILIY(func,Scalar,typename OtherDerived::Scalar);
00085 
00086   eigen_assert(size() == other.size());
00087 
00088   return internal::dot_nocheck<Derived,OtherDerived>::run(*this, other);
00089 }
00090 
00091 #ifdef EIGEN2_SUPPORT
00092 
00101 template<typename Derived>
00102 template<typename OtherDerived>
00103 typename internal::traits<Derived>::Scalar
00104 MatrixBase<Derived>::eigen2_dot(const MatrixBase<OtherDerived>& other) const
00105 {
00106   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
00107   EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
00108   EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
00109   EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value),
00110     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00111 
00112   eigen_assert(size() == other.size());
00113 
00114   return internal::dot_nocheck<OtherDerived,Derived>::run(other,*this);
00115 }
00116 #endif
00117 
00118 
00119 //---------- implementation of L2 norm and related functions ----------
00120 
00127 template<typename Derived>
00128 EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::squaredNorm() const
00129 {
00130   return internal::real((*this).cwiseAbs2().sum());
00131 }
00132 
00139 template<typename Derived>
00140 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::norm() const
00141 {
00142   return internal::sqrt(squaredNorm());
00143 }
00144 
00151 template<typename Derived>
00152 inline const typename MatrixBase<Derived>::PlainObject
00153 MatrixBase<Derived>::normalized() const
00154 {
00155   typedef typename internal::nested<Derived>::type Nested;
00156   typedef typename internal::remove_reference<Nested>::type _Nested;
00157   _Nested n(derived());
00158   return n / n.norm();
00159 }
00160 
00167 template<typename Derived>
00168 inline void MatrixBase<Derived>::normalize()
00169 {
00170   *this /= norm();
00171 }
00172 
00173 //---------- implementation of other norms ----------
00174 
00175 namespace internal {
00176 
00177 template<typename Derived, int p>
00178 struct lpNorm_selector
00179 {
00180   typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
00181   static inline RealScalar run(const MatrixBase<Derived>& m)
00182   {
00183     return pow(m.cwiseAbs().array().pow(p).sum(), RealScalar(1)/p);
00184   }
00185 };
00186 
00187 template<typename Derived>
00188 struct lpNorm_selector<Derived, 1>
00189 {
00190   static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
00191   {
00192     return m.cwiseAbs().sum();
00193   }
00194 };
00195 
00196 template<typename Derived>
00197 struct lpNorm_selector<Derived, 2>
00198 {
00199   static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
00200   {
00201     return m.norm();
00202   }
00203 };
00204 
00205 template<typename Derived>
00206 struct lpNorm_selector<Derived, Infinity>
00207 {
00208   static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(const MatrixBase<Derived>& m)
00209   {
00210     return m.cwiseAbs().maxCoeff();
00211   }
00212 };
00213 
00214 } // end namespace internal
00215 
00222 template<typename Derived>
00223 template<int p>
00224 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
00225 MatrixBase<Derived>::lpNorm() const
00226 {
00227   return internal::lpNorm_selector<Derived, p>::run(*this);
00228 }
00229 
00230 //---------- implementation of isOrthogonal / isUnitary ----------
00231 
00238 template<typename Derived>
00239 template<typename OtherDerived>
00240 bool MatrixBase<Derived>::isOrthogonal
00241 (const MatrixBase<OtherDerived>& other, RealScalar prec) const
00242 {
00243   typename internal::nested<Derived,2>::type nested(derived());
00244   typename internal::nested<OtherDerived,2>::type otherNested(other.derived());
00245   return internal::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm();
00246 }
00247 
00259 template<typename Derived>
00260 bool MatrixBase<Derived>::isUnitary(RealScalar prec) const
00261 {
00262   typename Derived::Nested nested(derived());
00263   for(Index i = 0; i < cols(); ++i)
00264   {
00265     if(!internal::isApprox(nested.col(i).squaredNorm(), static_cast<RealScalar>(1), prec))
00266       return false;
00267     for(Index j = 0; j < i; ++j)
00268       if(!internal::isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec))
00269         return false;
00270   }
00271   return true;
00272 }
00273 
00274 } // end namespace Eigen
00275 
00276 #endif // EIGEN_DOT_H