Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Core / Sum.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) 2008 Gael Guennebaud <g.gael@free.fr>
5 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
6 //
7 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25
26 #ifndef EIGEN_SUM_H
27 #define EIGEN_SUM_H
28
29 /***************************************************************************
30 * Part 1 : the logic deciding a strategy for vectorization and unrolling
31 ***************************************************************************/
32
33 template<typename Derived>
34 struct ei_sum_traits
35 {
36 private:
37   enum {
38     PacketSize = ei_packet_traits<typename Derived::Scalar>::size
39   };
40
41 public:
42   enum {
43     Vectorization = (int(Derived::Flags)&ActualPacketAccessBit)
44                  && (int(Derived::Flags)&LinearAccessBit)
45                   ? LinearVectorization
46                   : NoVectorization
47   };
48
49 private:
50   enum {
51     Cost = Derived::SizeAtCompileTime * Derived::CoeffReadCost
52            + (Derived::SizeAtCompileTime-1) * NumTraits<typename Derived::Scalar>::AddCost,
53     UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Vectorization) == int(NoVectorization) ? 1 : int(PacketSize))
54   };
55
56 public:
57   enum {
58     Unrolling = Cost <= UnrollingLimit
59               ? CompleteUnrolling
60               : NoUnrolling
61   };
62 };
63
64 /***************************************************************************
65 * Part 2 : unrollers
66 ***************************************************************************/
67
68 /*** no vectorization ***/
69
70 template<typename Derived, int Start, int Length>
71 struct ei_sum_novec_unroller
72 {
73   enum {
74     HalfLength = Length/2
75   };
76
77   typedef typename Derived::Scalar Scalar;
78
79   inline static Scalar run(const Derived &mat)
80   {
81     return ei_sum_novec_unroller<Derived, Start, HalfLength>::run(mat)
82          + ei_sum_novec_unroller<Derived, Start+HalfLength, Length-HalfLength>::run(mat);
83   }
84 };
85
86 template<typename Derived, int Start>
87 struct ei_sum_novec_unroller<Derived, Start, 1>
88 {
89   enum {
90     col = Start / Derived::RowsAtCompileTime,
91     row = Start % Derived::RowsAtCompileTime
92   };
93
94   typedef typename Derived::Scalar Scalar;
95
96   inline static Scalar run(const Derived &mat)
97   {
98     return mat.coeff(row, col);
99   }
100 };
101
102 /*** vectorization ***/
103   
104 template<typename Derived, int Start, int Length>
105 struct ei_sum_vec_unroller
106 {
107   enum {
108     PacketSize = ei_packet_traits<typename Derived::Scalar>::size,
109     HalfLength = Length/2
110   };
111
112   typedef typename Derived::Scalar Scalar;
113   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
114
115   inline static PacketScalar run(const Derived &mat)
116   {
117     return ei_padd(
118             ei_sum_vec_unroller<Derived, Start, HalfLength>::run(mat),
119             ei_sum_vec_unroller<Derived, Start+HalfLength, Length-HalfLength>::run(mat) );
120   }
121 };
122
123 template<typename Derived, int Start>
124 struct ei_sum_vec_unroller<Derived, Start, 1>
125 {
126   enum {
127     index = Start * ei_packet_traits<typename Derived::Scalar>::size,
128     row = int(Derived::Flags)&RowMajorBit
129         ? index / int(Derived::ColsAtCompileTime)
130         : index % Derived::RowsAtCompileTime,
131     col = int(Derived::Flags)&RowMajorBit
132         ? index % int(Derived::ColsAtCompileTime)
133         : index / Derived::RowsAtCompileTime,
134     alignment = (Derived::Flags & AlignedBit) ? Aligned : Unaligned
135   };
136
137   typedef typename Derived::Scalar Scalar;
138   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
139
140   inline static PacketScalar run(const Derived &mat)
141   {
142     return mat.template packet<alignment>(row, col);
143   }
144 };
145
146 /***************************************************************************
147 * Part 3 : implementation of all cases
148 ***************************************************************************/
149
150 template<typename Derived,
151          int Vectorization = ei_sum_traits<Derived>::Vectorization,
152          int Unrolling = ei_sum_traits<Derived>::Unrolling
153 >
154 struct ei_sum_impl;
155
156 template<typename Derived>
157 struct ei_sum_impl<Derived, NoVectorization, NoUnrolling>
158 {
159   typedef typename Derived::Scalar Scalar;
160   static Scalar run(const Derived& mat)
161   {
162     ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using a non initialized matrix");
163     Scalar res;
164     res = mat.coeff(0, 0);
165     for(int i = 1; i < mat.rows(); ++i)
166       res += mat.coeff(i, 0);
167     for(int j = 1; j < mat.cols(); ++j)
168       for(int i = 0; i < mat.rows(); ++i)
169         res += mat.coeff(i, j);
170     return res;
171   }
172 };
173
174 template<typename Derived>
175 struct ei_sum_impl<Derived, NoVectorization, CompleteUnrolling>
176   : public ei_sum_novec_unroller<Derived, 0, Derived::SizeAtCompileTime>
177 {};
178
179 template<typename Derived>
180 struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling>
181 {
182   typedef typename Derived::Scalar Scalar;
183   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
184
185   static Scalar run(const Derived& mat)
186   {
187     const int size = mat.size();
188     const int packetSize = ei_packet_traits<Scalar>::size;
189     const int alignedStart =  (Derived::Flags & AlignedBit)
190                            || !(Derived::Flags & DirectAccessBit)
191                            ? 0
192                            : ei_alignmentOffset(&mat.const_cast_derived().coeffRef(0), size);
193     enum {
194       alignment = (Derived::Flags & DirectAccessBit) || (Derived::Flags & AlignedBit)
195                 ? Aligned : Unaligned
196     };
197     const int alignedSize = ((size-alignedStart)/packetSize)*packetSize;
198     const int alignedEnd = alignedStart + alignedSize;
199     Scalar res;
200
201     if(alignedSize)
202     {
203       PacketScalar packet_res = mat.template packet<alignment>(alignedStart);
204       for(int index = alignedStart + packetSize; index < alignedEnd; index += packetSize)
205         packet_res = ei_padd(packet_res, mat.template packet<alignment>(index));
206       res = ei_predux(packet_res);
207     }
208     else // too small to vectorize anything.
209          // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
210     {
211       res = Scalar(0);
212     }
213
214     for(int index = 0; index < alignedStart; ++index)
215       res += mat.coeff(index);
216
217     for(int index = alignedEnd; index < size; ++index)
218       res += mat.coeff(index);
219
220     return res;
221   }
222 };
223
224 template<typename Derived>
225 struct ei_sum_impl<Derived, LinearVectorization, CompleteUnrolling>
226 {
227   typedef typename Derived::Scalar Scalar;
228   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
229   enum {
230     PacketSize = ei_packet_traits<Scalar>::size,
231     Size = Derived::SizeAtCompileTime,
232     VectorizationSize = (Size / PacketSize) * PacketSize
233   };
234   static Scalar run(const Derived& mat)
235   {
236     Scalar res = ei_predux(ei_sum_vec_unroller<Derived, 0, Size / PacketSize>::run(mat));
237     if (VectorizationSize != Size)
238       res += ei_sum_novec_unroller<Derived, VectorizationSize, Size-VectorizationSize>::run(mat);
239     return res;
240   }
241 };
242
243 /***************************************************************************
244 * Part 4 : implementation of MatrixBase methods
245 ***************************************************************************/
246
247 /** \returns the sum of all coefficients of *this
248   *
249   * \sa trace()
250   */
251 template<typename Derived>
252 inline typename ei_traits<Derived>::Scalar
253 MatrixBase<Derived>::sum() const
254 {
255   return ei_sum_impl<Derived>::run(derived());
256 }
257
258 /** \returns the trace of \c *this, i.e. the sum of the coefficients on the main diagonal.
259   *
260   * \c *this can be any matrix, not necessarily square.
261   *
262   * \sa diagonal(), sum()
263   */
264 template<typename Derived>
265 inline typename ei_traits<Derived>::Scalar
266 MatrixBase<Derived>::trace() const
267 {
268   return diagonal().sum();
269 }
270
271 #endif // EIGEN_SUM_H