Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / extern / Eigen2 / Eigen / src / Core / Dot.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra. Eigen itself is part of the KDE project.
3 //
4 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24
25 #ifndef EIGEN_DOT_H
26 #define EIGEN_DOT_H
27
28 /***************************************************************************
29 * Part 1 : the logic deciding a strategy for vectorization and unrolling
30 ***************************************************************************/
31
32 template<typename Derived1, typename Derived2>
33 struct ei_dot_traits
34 {
35 public:
36   enum {
37     Vectorization = (int(Derived1::Flags)&int(Derived2::Flags)&ActualPacketAccessBit)
38                  && (int(Derived1::Flags)&int(Derived2::Flags)&LinearAccessBit)
39                   ? LinearVectorization
40                   : NoVectorization
41   };
42
43 private:
44   typedef typename Derived1::Scalar Scalar;
45   enum {
46     PacketSize = ei_packet_traits<Scalar>::size,
47     Cost = Derived1::SizeAtCompileTime * (Derived1::CoeffReadCost + Derived2::CoeffReadCost + NumTraits<Scalar>::MulCost)
48            + (Derived1::SizeAtCompileTime-1) * NumTraits<Scalar>::AddCost,
49     UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Vectorization) == int(NoVectorization) ? 1 : int(PacketSize))
50   };
51
52 public:
53   enum {
54     Unrolling = Cost <= UnrollingLimit
55               ? CompleteUnrolling
56               : NoUnrolling
57   };
58 };
59
60 /***************************************************************************
61 * Part 2 : unrollers
62 ***************************************************************************/
63
64 /*** no vectorization ***/
65
66 template<typename Derived1, typename Derived2, int Start, int Length>
67 struct ei_dot_novec_unroller
68 {
69   enum {
70     HalfLength = Length/2
71   };
72
73   typedef typename Derived1::Scalar Scalar;
74
75   inline static Scalar run(const Derived1& v1, const Derived2& v2)
76   {
77     return ei_dot_novec_unroller<Derived1, Derived2, Start, HalfLength>::run(v1, v2)
78          + ei_dot_novec_unroller<Derived1, Derived2, Start+HalfLength, Length-HalfLength>::run(v1, v2);
79   }
80 };
81
82 template<typename Derived1, typename Derived2, int Start>
83 struct ei_dot_novec_unroller<Derived1, Derived2, Start, 1>
84 {
85   typedef typename Derived1::Scalar Scalar;
86
87   inline static Scalar run(const Derived1& v1, const Derived2& v2)
88   {
89     return v1.coeff(Start) * ei_conj(v2.coeff(Start));
90   }
91 };
92
93 /*** vectorization ***/
94
95 template<typename Derived1, typename Derived2, int Index, int Stop,
96          bool LastPacket = (Stop-Index == ei_packet_traits<typename Derived1::Scalar>::size)>
97 struct ei_dot_vec_unroller
98 {
99   typedef typename Derived1::Scalar Scalar;
100   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
101
102   enum {
103     row1 = Derived1::RowsAtCompileTime == 1 ? 0 : Index,
104     col1 = Derived1::RowsAtCompileTime == 1 ? Index : 0,
105     row2 = Derived2::RowsAtCompileTime == 1 ? 0 : Index,
106     col2 = Derived2::RowsAtCompileTime == 1 ? Index : 0
107   };
108
109   inline static PacketScalar run(const Derived1& v1, const Derived2& v2)
110   {
111     return ei_pmadd(
112       v1.template packet<Aligned>(row1, col1),
113       v2.template packet<Aligned>(row2, col2),
114       ei_dot_vec_unroller<Derived1, Derived2, Index+ei_packet_traits<Scalar>::size, Stop>::run(v1, v2)
115     );
116   }
117 };
118
119 template<typename Derived1, typename Derived2, int Index, int Stop>
120 struct ei_dot_vec_unroller<Derived1, Derived2, Index, Stop, true>
121 {
122   enum {
123     row1 = Derived1::RowsAtCompileTime == 1 ? 0 : Index,
124     col1 = Derived1::RowsAtCompileTime == 1 ? Index : 0,
125     row2 = Derived2::RowsAtCompileTime == 1 ? 0 : Index,
126     col2 = Derived2::RowsAtCompileTime == 1 ? Index : 0,
127     alignment1 = (Derived1::Flags & AlignedBit) ? Aligned : Unaligned,
128     alignment2 = (Derived2::Flags & AlignedBit) ? Aligned : Unaligned
129   };
130
131   typedef typename Derived1::Scalar Scalar;
132   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
133
134   inline static PacketScalar run(const Derived1& v1, const Derived2& v2)
135   {
136     return ei_pmul(v1.template packet<alignment1>(row1, col1), v2.template packet<alignment2>(row2, col2));
137   }
138 };
139
140 /***************************************************************************
141 * Part 3 : implementation of all cases
142 ***************************************************************************/
143
144 template<typename Derived1, typename Derived2,
145          int Vectorization = ei_dot_traits<Derived1, Derived2>::Vectorization,
146          int Unrolling = ei_dot_traits<Derived1, Derived2>::Unrolling
147 >
148 struct ei_dot_impl;
149
150 template<typename Derived1, typename Derived2>
151 struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling>
152 {
153   typedef typename Derived1::Scalar Scalar;
154   static Scalar run(const Derived1& v1, const Derived2& v2)
155   {
156     ei_assert(v1.size()>0 && "you are using a non initialized vector");
157     Scalar res;
158     res = v1.coeff(0) * ei_conj(v2.coeff(0));
159     for(int i = 1; i < v1.size(); ++i)
160       res += v1.coeff(i) * ei_conj(v2.coeff(i));
161     return res;
162   }
163 };
164
165 template<typename Derived1, typename Derived2>
166 struct ei_dot_impl<Derived1, Derived2, NoVectorization, CompleteUnrolling>
167   : public ei_dot_novec_unroller<Derived1, Derived2, 0, Derived1::SizeAtCompileTime>
168 {};
169
170 template<typename Derived1, typename Derived2>
171 struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>
172 {
173   typedef typename Derived1::Scalar Scalar;
174   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
175
176   static Scalar run(const Derived1& v1, const Derived2& v2)
177   {
178     const int size = v1.size();
179     const int packetSize = ei_packet_traits<Scalar>::size;
180     const int alignedSize = (size/packetSize)*packetSize;
181     enum {
182       alignment1 = (Derived1::Flags & AlignedBit) ? Aligned : Unaligned,
183       alignment2 = (Derived2::Flags & AlignedBit) ? Aligned : Unaligned
184     };
185     Scalar res;
186
187     // do the vectorizable part of the sum
188     if(size >= packetSize)
189     {
190       PacketScalar packet_res = ei_pmul(
191                                   v1.template packet<alignment1>(0),
192                                   v2.template packet<alignment2>(0)
193                                 );
194       for(int index = packetSize; index<alignedSize; index += packetSize)
195       {
196         packet_res = ei_pmadd(
197                        v1.template packet<alignment1>(index),
198                        v2.template packet<alignment2>(index),
199                        packet_res
200                      );
201       }
202       res = ei_predux(packet_res);
203
204       // now we must do the rest without vectorization.
205       if(alignedSize == size) return res;
206     }
207     else // too small to vectorize anything.
208          // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
209     {
210       res = Scalar(0);
211     }
212
213     // do the remainder of the vector
214     for(int index = alignedSize; index < size; ++index)
215     {
216       res += v1.coeff(index) * v2.coeff(index);
217     }
218
219     return res;
220   }
221 };
222
223 template<typename Derived1, typename Derived2>
224 struct ei_dot_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling>
225 {
226   typedef typename Derived1::Scalar Scalar;
227   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
228   enum {
229     PacketSize = ei_packet_traits<Scalar>::size,
230     Size = Derived1::SizeAtCompileTime,
231     VectorizationSize = (Size / PacketSize) * PacketSize
232   };
233   static Scalar run(const Derived1& v1, const Derived2& v2)
234   {
235     Scalar res = ei_predux(ei_dot_vec_unroller<Derived1, Derived2, 0, VectorizationSize>::run(v1, v2));
236     if (VectorizationSize != Size)
237       res += ei_dot_novec_unroller<Derived1, Derived2, VectorizationSize, Size-VectorizationSize>::run(v1, v2);
238     return res;
239   }
240 };
241
242 /***************************************************************************
243 * Part 4 : implementation of MatrixBase methods
244 ***************************************************************************/
245
246 /** \returns the dot product of *this with other.
247   *
248   * \only_for_vectors
249   *
250   * \note If the scalar type is complex numbers, then this function returns the hermitian
251   * (sesquilinear) dot product, linear in the first variable and conjugate-linear in the
252   * second variable.
253   *
254   * \sa squaredNorm(), norm()
255   */
256 template<typename Derived>
257 template<typename OtherDerived>
258 typename ei_traits<Derived>::Scalar
259 MatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const
260 {
261   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
262   EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
263   EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
264   EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret),
265     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
266
267   ei_assert(size() == other.size());
268
269   return ei_dot_impl<Derived, OtherDerived>::run(derived(), other.derived());
270 }
271
272 /** \returns the squared \em l2 norm of *this, i.e., for vectors, the dot product of *this with itself.
273   *
274   * \sa dot(), norm()
275   */
276 template<typename Derived>
277 inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<Derived>::squaredNorm() const
278 {
279   return ei_real((*this).cwise().abs2().sum());
280 }
281
282 /** \returns the \em l2 norm of *this, i.e., for vectors, the square root of the dot product of *this with itself.
283   *
284   * \sa dot(), squaredNorm()
285   */
286 template<typename Derived>
287 inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<Derived>::norm() const
288 {
289   return ei_sqrt(squaredNorm());
290 }
291
292 /** \returns an expression of the quotient of *this by its own norm.
293   *
294   * \only_for_vectors
295   *
296   * \sa norm(), normalize()
297   */
298 template<typename Derived>
299 inline const typename MatrixBase<Derived>::PlainMatrixType
300 MatrixBase<Derived>::normalized() const
301 {
302   typedef typename ei_nested<Derived>::type Nested;
303   typedef typename ei_unref<Nested>::type _Nested;
304   _Nested n(derived());
305   return n / n.norm();
306 }
307
308 /** Normalizes the vector, i.e. divides it by its own norm.
309   *
310   * \only_for_vectors
311   *
312   * \sa norm(), normalized()
313   */
314 template<typename Derived>
315 inline void MatrixBase<Derived>::normalize()
316 {
317   *this /= norm();
318 }
319
320 /** \returns true if *this is approximately orthogonal to \a other,
321   *          within the precision given by \a prec.
322   *
323   * Example: \include MatrixBase_isOrthogonal.cpp
324   * Output: \verbinclude MatrixBase_isOrthogonal.out
325   */
326 template<typename Derived>
327 template<typename OtherDerived>
328 bool MatrixBase<Derived>::isOrthogonal
329 (const MatrixBase<OtherDerived>& other, RealScalar prec) const
330 {
331   typename ei_nested<Derived,2>::type nested(derived());
332   typename ei_nested<OtherDerived,2>::type otherNested(other.derived());
333   return ei_abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm();
334 }
335
336 /** \returns true if *this is approximately an unitary matrix,
337   *          within the precision given by \a prec. In the case where the \a Scalar
338   *          type is real numbers, a unitary matrix is an orthogonal matrix, whence the name.
339   *
340   * \note This can be used to check whether a family of vectors forms an orthonormal basis.
341   *       Indeed, \c m.isUnitary() returns true if and only if the columns (equivalently, the rows) of m form an
342   *       orthonormal basis.
343   *
344   * Example: \include MatrixBase_isUnitary.cpp
345   * Output: \verbinclude MatrixBase_isUnitary.out
346   */
347 template<typename Derived>
348 bool MatrixBase<Derived>::isUnitary(RealScalar prec) const
349 {
350   typename Derived::Nested nested(derived());
351   for(int i = 0; i < cols(); ++i)
352   {
353     if(!ei_isApprox(nested.col(i).squaredNorm(), static_cast<Scalar>(1), prec))
354       return false;
355     for(int j = 0; j < i; ++j)
356       if(!ei_isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec))
357         return false;
358   }
359   return true;
360 }
361 #endif // EIGEN_DOT_H