vcglib: Dot.h Source File

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2006-2008 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 /***************************************************************************
00029 * Part 1 : the logic deciding a strategy for vectorization and unrolling
00030 ***************************************************************************/
00031 
00032 template<typename Derived1, typename Derived2>
00033 struct ei_dot_traits
00034 {
00035 public:
00036   enum {
00037     Vectorization = (int(Derived1::Flags)&int(Derived2::Flags)&ActualPacketAccessBit)
00038                  && (int(Derived1::Flags)&int(Derived2::Flags)&LinearAccessBit)
00039                   ? LinearVectorization
00040                   : NoVectorization
00041   };
00042 
00043 private:
00044   typedef typename Derived1::Scalar Scalar;
00045   enum {
00046     PacketSize = ei_packet_traits<Scalar>::size,
00047     Cost = Derived1::SizeAtCompileTime * (Derived1::CoeffReadCost + Derived2::CoeffReadCost + NumTraits<Scalar>::MulCost)
00048            + (Derived1::SizeAtCompileTime-1) * NumTraits<Scalar>::AddCost,
00049     UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Vectorization) == int(NoVectorization) ? 1 : int(PacketSize))
00050   };
00051 
00052 public:
00053   enum {
00054     Unrolling = Cost <= UnrollingLimit
00055               ? CompleteUnrolling
00056               : NoUnrolling
00057   };
00058 };
00059 
00060 /***************************************************************************
00061 * Part 2 : unrollers
00062 ***************************************************************************/
00063 
00064 /*** no vectorization ***/
00065 
00066 template<typename Derived1, typename Derived2, int Start, int Length>
00067 struct ei_dot_novec_unroller
00068 {
00069   enum {
00070     HalfLength = Length/2
00071   };
00072 
00073   typedef typename Derived1::Scalar Scalar;
00074 
00075   inline static Scalar run(const Derived1& v1, const Derived2& v2)
00076   {
00077     return ei_dot_novec_unroller<Derived1, Derived2, Start, HalfLength>::run(v1, v2)
00078          + ei_dot_novec_unroller<Derived1, Derived2, Start+HalfLength, Length-HalfLength>::run(v1, v2);
00079   }
00080 };
00081 
00082 template<typename Derived1, typename Derived2, int Start>
00083 struct ei_dot_novec_unroller<Derived1, Derived2, Start, 1>
00084 {
00085   typedef typename Derived1::Scalar Scalar;
00086 
00087   inline static Scalar run(const Derived1& v1, const Derived2& v2)
00088   {
00089     return v1.coeff(Start) * ei_conj(v2.coeff(Start));
00090   }
00091 };
00092 
00093 /*** vectorization ***/
00094 
00095 template<typename Derived1, typename Derived2, int Index, int Stop,
00096          bool LastPacket = (Stop-Index == ei_packet_traits<typename Derived1::Scalar>::size)>
00097 struct ei_dot_vec_unroller
00098 {
00099   typedef typename Derived1::Scalar Scalar;
00100   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
00101 
00102   enum {
00103     row1 = Derived1::RowsAtCompileTime == 1 ? 0 : Index,
00104     col1 = Derived1::RowsAtCompileTime == 1 ? Index : 0,
00105     row2 = Derived2::RowsAtCompileTime == 1 ? 0 : Index,
00106     col2 = Derived2::RowsAtCompileTime == 1 ? Index : 0
00107   };
00108 
00109   inline static PacketScalar run(const Derived1& v1, const Derived2& v2)
00110   {
00111     return ei_pmadd(
00112       v1.template packet<Aligned>(row1, col1),
00113       v2.template packet<Aligned>(row2, col2),
00114       ei_dot_vec_unroller<Derived1, Derived2, Index+ei_packet_traits<Scalar>::size, Stop>::run(v1, v2)
00115     );
00116   }
00117 };
00118 
00119 template<typename Derived1, typename Derived2, int Index, int Stop>
00120 struct ei_dot_vec_unroller<Derived1, Derived2, Index, Stop, true>
00121 {
00122   enum {
00123     row1 = Derived1::RowsAtCompileTime == 1 ? 0 : Index,
00124     col1 = Derived1::RowsAtCompileTime == 1 ? Index : 0,
00125     row2 = Derived2::RowsAtCompileTime == 1 ? 0 : Index,
00126     col2 = Derived2::RowsAtCompileTime == 1 ? Index : 0,
00127     alignment1 = (Derived1::Flags & AlignedBit) ? Aligned : Unaligned,
00128     alignment2 = (Derived2::Flags & AlignedBit) ? Aligned : Unaligned
00129   };
00130 
00131   typedef typename Derived1::Scalar Scalar;
00132   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
00133 
00134   inline static PacketScalar run(const Derived1& v1, const Derived2& v2)
00135   {
00136     return ei_pmul(v1.template packet<alignment1>(row1, col1), v2.template packet<alignment2>(row2, col2));
00137   }
00138 };
00139 
00140 /***************************************************************************
00141 * Part 3 : implementation of all cases
00142 ***************************************************************************/
00143 
00144 template<typename Derived1, typename Derived2,
00145          int Vectorization = ei_dot_traits<Derived1, Derived2>::Vectorization,
00146          int Unrolling = ei_dot_traits<Derived1, Derived2>::Unrolling
00147 >
00148 struct ei_dot_impl;
00149 
00150 template<typename Derived1, typename Derived2>
00151 struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling>
00152 {
00153   typedef typename Derived1::Scalar Scalar;
00154   static Scalar run(const Derived1& v1, const Derived2& v2)
00155   {
00156     ei_assert(v1.size()>0 && "you are using a non initialized vector");
00157     Scalar res;
00158     res = v1.coeff(0) * ei_conj(v2.coeff(0));
00159     for(int i = 1; i < v1.size(); ++i)
00160       res += v1.coeff(i) * ei_conj(v2.coeff(i));
00161     return res;
00162   }
00163 };
00164 
00165 template<typename Derived1, typename Derived2>
00166 struct ei_dot_impl<Derived1, Derived2, NoVectorization, CompleteUnrolling>
00167   : public ei_dot_novec_unroller<Derived1, Derived2, 0, Derived1::SizeAtCompileTime>
00168 {};
00169 
00170 template<typename Derived1, typename Derived2>
00171 struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>
00172 {
00173   typedef typename Derived1::Scalar Scalar;
00174   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
00175 
00176   static Scalar run(const Derived1& v1, const Derived2& v2)
00177   {
00178     const int size = v1.size();
00179     const int packetSize = ei_packet_traits<Scalar>::size;
00180     const int alignedSize = (size/packetSize)*packetSize;
00181     enum {
00182       alignment1 = (Derived1::Flags & AlignedBit) ? Aligned : Unaligned,
00183       alignment2 = (Derived2::Flags & AlignedBit) ? Aligned : Unaligned
00184     };
00185     Scalar res;
00186 
00187     // do the vectorizable part of the sum
00188     if(size >= packetSize)
00189     {
00190       PacketScalar packet_res = ei_pmul(
00191                                   v1.template packet<alignment1>(0),
00192                                   v2.template packet<alignment2>(0)
00193                                 );
00194       for(int index = packetSize; index<alignedSize; index += packetSize)
00195       {
00196         packet_res = ei_pmadd(
00197                        v1.template packet<alignment1>(index),
00198                        v2.template packet<alignment2>(index),
00199                        packet_res
00200                      );
00201       }
00202       res = ei_predux(packet_res);
00203 
00204       // now we must do the rest without vectorization.
00205       if(alignedSize == size) return res;
00206     }
00207     else // too small to vectorize anything.
00208          // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
00209     {
00210       res = Scalar(0);
00211     }
00212 
00213     // do the remainder of the vector
00214     for(int index = alignedSize; index < size; ++index)
00215     {
00216       res += v1.coeff(index) * v2.coeff(index);
00217     }
00218 
00219     return res;
00220   }
00221 };
00222 
00223 template<typename Derived1, typename Derived2>
00224 struct ei_dot_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling>
00225 {
00226   typedef typename Derived1::Scalar Scalar;
00227   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
00228   enum {
00229     PacketSize = ei_packet_traits<Scalar>::size,
00230     Size = Derived1::SizeAtCompileTime,
00231     VectorizationSize = (Size / PacketSize) * PacketSize
00232   };
00233   static Scalar run(const Derived1& v1, const Derived2& v2)
00234   {
00235     Scalar res = ei_predux(ei_dot_vec_unroller<Derived1, Derived2, 0, VectorizationSize>::run(v1, v2));
00236     if (VectorizationSize != Size)
00237       res += ei_dot_novec_unroller<Derived1, Derived2, VectorizationSize, Size-VectorizationSize>::run(v1, v2);
00238     return res;
00239   }
00240 };
00241 
00242 /***************************************************************************
00243 * Part 4 : implementation of MatrixBase methods
00244 ***************************************************************************/
00245 
00256 template<typename Derived>
00257 template<typename OtherDerived>
00258 typename ei_traits<Derived>::Scalar
00259 MatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const
00260 {
00261   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
00262   EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
00263   EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
00264   EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret),
00265     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00266 
00267   ei_assert(size() == other.size());
00268 
00269   return ei_dot_impl<Derived, OtherDerived>::run(derived(), other.derived());
00270 }
00271 
00276 template<typename Derived>
00277 inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<Derived>::squaredNorm() const
00278 {
00279   return ei_real((*this).cwise().abs2().sum());
00280 }
00281 
00286 template<typename Derived>
00287 inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<Derived>::norm() const
00288 {
00289   return ei_sqrt(squaredNorm());
00290 }
00291 
00298 template<typename Derived>
00299 inline const typename MatrixBase<Derived>::PlainMatrixType
00300 MatrixBase<Derived>::normalized() const
00301 {
00302   typedef typename ei_nested<Derived>::type Nested;
00303   typedef typename ei_unref<Nested>::type _Nested;
00304   _Nested n(derived());
00305   return n / n.norm();
00306 }
00307 
00314 template<typename Derived>
00315 inline void MatrixBase<Derived>::normalize()
00316 {
00317   *this /= norm();
00318 }
00319 
00326 template<typename Derived>
00327 template<typename OtherDerived>
00328 bool MatrixBase<Derived>::isOrthogonal
00329 (const MatrixBase<OtherDerived>& other, RealScalar prec) const
00330 {
00331   typename ei_nested<Derived,2>::type nested(derived());
00332   typename ei_nested<OtherDerived,2>::type otherNested(other.derived());
00333   return ei_abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm();
00334 }
00335 
00347 template<typename Derived>
00348 bool MatrixBase<Derived>::isUnitary(RealScalar prec) const
00349 {
00350   typename Derived::Nested nested(derived());
00351   for(int i = 0; i < cols(); ++i)
00352   {
00353     if(!ei_isApprox(nested.col(i).squaredNorm(), static_cast<Scalar>(1), prec))
00354       return false;
00355     for(int j = 0; j < i; ++j)
00356       if(!ei_isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec))
00357         return false;
00358   }
00359   return true;
00360 }
00361 #endif // EIGEN_DOT_H