Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Core / arch / AltiVec / PacketMath.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 Konstantinos Margaritis <markos@codex.gr>
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_PACKET_MATH_ALTIVEC_H
26 #define EIGEN_PACKET_MATH_ALTIVEC_H
27
28 #ifndef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
29 #define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 4
30 #endif
31
32 typedef __vector float          v4f;
33 typedef __vector int            v4i;
34 typedef __vector unsigned int   v4ui;
35 typedef __vector __bool int     v4bi;
36
37 // We don't want to write the same code all the time, but we need to reuse the constants
38 // and it doesn't really work to declare them global, so we define macros instead
39
40 #define USE_CONST_v0i     const v4i   v0i   = vec_splat_s32(0)
41 #define USE_CONST_v1i     const v4i   v1i   = vec_splat_s32(1)
42 #define USE_CONST_v16i_   const v4i   v16i_ = vec_splat_s32(-16)
43 #define USE_CONST_v0f     USE_CONST_v0i; const v4f v0f = (v4f) v0i
44 #define USE_CONST_v1f     USE_CONST_v1i; const v4f v1f = vec_ctf(v1i, 0)
45 #define USE_CONST_v1i_    const v4ui  v1i_  = vec_splat_u32(-1)
46 #define USE_CONST_v0f_    USE_CONST_v1i_; const v4f v0f_ = (v4f) vec_sl(v1i_, v1i_)
47
48 template<> struct ei_packet_traits<float>  { typedef v4f type; enum {size=4}; };
49 template<> struct ei_packet_traits<int>    { typedef v4i type; enum {size=4}; };
50
51 template<> struct ei_unpacket_traits<v4f>  { typedef float  type; enum {size=4}; };
52 template<> struct ei_unpacket_traits<v4i>  { typedef int    type; enum {size=4}; };
53
54 inline std::ostream & operator <<(std::ostream & s, const v4f & v)
55 {
56   union {
57     v4f   v;
58     float n[4];
59   } vt;
60   vt.v = v;
61   s << vt.n[0] << ", " << vt.n[1] << ", " << vt.n[2] << ", " << vt.n[3];
62   return s;
63 }
64
65 inline std::ostream & operator <<(std::ostream & s, const v4i & v)
66 {
67   union {
68     v4i   v;
69     int n[4];
70   } vt;
71   vt.v = v;
72   s << vt.n[0] << ", " << vt.n[1] << ", " << vt.n[2] << ", " << vt.n[3];
73   return s;
74 }
75
76 inline std::ostream & operator <<(std::ostream & s, const v4ui & v)
77 {
78   union {
79     v4ui   v;
80     unsigned int n[4];
81   } vt;
82   vt.v = v;
83   s << vt.n[0] << ", " << vt.n[1] << ", " << vt.n[2] << ", " << vt.n[3];
84   return s;
85 }
86
87 inline std::ostream & operator <<(std::ostream & s, const v4bi & v)
88 {
89   union {
90     __vector __bool int v;
91     unsigned int n[4];
92   } vt;
93   vt.v = v;
94   s << vt.n[0] << ", " << vt.n[1] << ", " << vt.n[2] << ", " << vt.n[3];
95   return s;
96 }
97
98 template<> inline v4f  ei_padd(const v4f&   a, const v4f&   b) { return vec_add(a,b); }
99 template<> inline v4i  ei_padd(const v4i&   a, const v4i&   b) { return vec_add(a,b); }
100
101 template<> inline v4f  ei_psub(const v4f&   a, const v4f&   b) { return vec_sub(a,b); }
102 template<> inline v4i  ei_psub(const v4i&   a, const v4i&   b) { return vec_sub(a,b); }
103
104 template<> inline v4f  ei_pmul(const v4f&   a, const v4f&   b) { USE_CONST_v0f; return vec_madd(a,b, v0f); }
105 template<> inline v4i  ei_pmul(const v4i&   a, const v4i&   b)
106 {
107   // Detailed in: http://freevec.org/content/32bit_signed_integer_multiplication_altivec
108   //Set up constants, variables
109   v4i a1, b1, bswap, low_prod, high_prod, prod, prod_, v1sel;
110   USE_CONST_v0i;
111   USE_CONST_v1i;
112   USE_CONST_v16i_;
113
114   // Get the absolute values 
115   a1  = vec_abs(a);
116   b1  = vec_abs(b);
117
118   // Get the signs using xor
119   v4bi sgn = (v4bi) vec_cmplt(vec_xor(a, b), v0i);
120
121   // Do the multiplication for the asbolute values.
122   bswap = (v4i) vec_rl((v4ui) b1, (v4ui) v16i_ );
123   low_prod = vec_mulo((__vector short)a1, (__vector short)b1);
124   high_prod = vec_msum((__vector short)a1, (__vector short)bswap, v0i);
125   high_prod = (v4i) vec_sl((v4ui) high_prod, (v4ui) v16i_);
126   prod = vec_add( low_prod, high_prod );
127
128   // NOR the product and select only the negative elements according to the sign mask
129   prod_ = vec_nor(prod, prod);
130   prod_ = vec_sel(v0i, prod_, sgn);
131
132   // Add 1 to the result to get the negative numbers
133   v1sel = vec_sel(v0i, v1i, sgn);
134   prod_ = vec_add(prod_, v1sel);
135
136   // Merge the results back to the final vector.
137   prod = vec_sel(prod, prod_, sgn);
138
139   return prod;
140 }
141
142 template<> inline v4f  ei_pdiv(const v4f&   a, const v4f&   b) {
143   v4f t, y_0, y_1, res;
144   USE_CONST_v0f;
145   USE_CONST_v1f;
146
147   // Altivec does not offer a divide instruction, we have to do a reciprocal approximation
148   y_0 = vec_re(b);
149   
150   // Do one Newton-Raphson iteration to get the needed accuracy
151   t = vec_nmsub(y_0, b, v1f);
152   y_1 = vec_madd(y_0, t, y_0);
153
154   res = vec_madd(a, y_1, v0f);
155   return res;
156 }
157
158 template<> inline v4f  ei_pmadd(const v4f&  a, const v4f&   b, const v4f&  c) { return vec_madd(a, b, c); }
159
160 template<> inline v4f  ei_pmin(const v4f&   a, const v4f&   b) { return vec_min(a,b); }
161 template<> inline v4i  ei_pmin(const v4i&   a, const v4i&   b) { return vec_min(a,b); }
162
163 template<> inline v4f  ei_pmax(const v4f&   a, const v4f&   b) { return vec_max(a,b); }
164 template<> inline v4i  ei_pmax(const v4i&   a, const v4i&   b) { return vec_max(a,b); }
165
166 template<> inline v4f  ei_pload(const float* from) { return vec_ld(0, from); }
167 template<> inline v4i  ei_pload(const int*   from) { return vec_ld(0, from); }
168
169 template<> inline v4f  ei_ploadu(const float*  from)
170 {
171   // Taken from http://developer.apple.com/hardwaredrivers/ve/alignment.html
172   __vector unsigned char MSQ, LSQ;
173   __vector unsigned char mask;
174   MSQ = vec_ld(0, (unsigned char *)from);          // most significant quadword
175   LSQ = vec_ld(15, (unsigned char *)from);         // least significant quadword
176   mask = vec_lvsl(0, from);                        // create the permute mask
177   return (v4f) vec_perm(MSQ, LSQ, mask);           // align the data
178 }
179
180 template<> inline v4i    ei_ploadu(const int*    from)
181 {
182   // Taken from http://developer.apple.com/hardwaredrivers/ve/alignment.html
183   __vector unsigned char MSQ, LSQ;
184   __vector unsigned char mask;
185   MSQ = vec_ld(0, (unsigned char *)from);          // most significant quadword
186   LSQ = vec_ld(15, (unsigned char *)from);         // least significant quadword
187   mask = vec_lvsl(0, from);                        // create the permute mask
188   return (v4i) vec_perm(MSQ, LSQ, mask);    // align the data
189 }
190
191 template<> inline v4f  ei_pset1(const float&  from)
192 {
193   // Taken from http://developer.apple.com/hardwaredrivers/ve/alignment.html
194   float __attribute__(aligned(16)) af[4];
195   af[0] = from;
196   v4f vc = vec_ld(0, af);
197   vc = vec_splat(vc, 0);
198   return vc;
199 }
200
201 template<> inline v4i    ei_pset1(const int&    from)
202 {
203   int __attribute__(aligned(16)) ai[4];
204   ai[0] = from;
205   v4i vc = vec_ld(0, ai);
206   vc = vec_splat(vc, 0);
207   return vc;
208 }
209
210 template<> inline void ei_pstore(float*   to, const v4f&   from) { vec_st(from, 0, to); }
211 template<> inline void ei_pstore(int*     to, const v4i&   from) { vec_st(from, 0, to); }
212
213 template<> inline void ei_pstoreu(float*  to, const v4f&   from)
214 {
215   // Taken from http://developer.apple.com/hardwaredrivers/ve/alignment.html
216   // Warning: not thread safe!
217   __vector unsigned char MSQ, LSQ, edges;
218   __vector unsigned char edgeAlign, align;
219
220   MSQ = vec_ld(0, (unsigned char *)to);                     // most significant quadword
221   LSQ = vec_ld(15, (unsigned char *)to);                    // least significant quadword
222   edgeAlign = vec_lvsl(0, to);                              // permute map to extract edges
223   edges=vec_perm(LSQ,MSQ,edgeAlign);                        // extract the edges
224   align = vec_lvsr( 0, to );                                // permute map to misalign data
225   MSQ = vec_perm(edges,(__vector unsigned char)from,align);   // misalign the data (MSQ)
226   LSQ = vec_perm((__vector unsigned char)from,edges,align);   // misalign the data (LSQ)
227   vec_st( LSQ, 15, (unsigned char *)to );                   // Store the LSQ part first
228   vec_st( MSQ, 0, (unsigned char *)to );                    // Store the MSQ part
229 }
230
231 template<> inline void ei_pstoreu(int*    to , const v4i&    from )
232 {
233   // Taken from http://developer.apple.com/hardwaredrivers/ve/alignment.html
234   // Warning: not thread safe!
235   __vector unsigned char MSQ, LSQ, edges;
236   __vector unsigned char edgeAlign, align;
237
238   MSQ = vec_ld(0, (unsigned char *)to);                     // most significant quadword
239   LSQ = vec_ld(15, (unsigned char *)to);                    // least significant quadword
240   edgeAlign = vec_lvsl(0, to);                              // permute map to extract edges
241   edges=vec_perm(LSQ,MSQ,edgeAlign);                        // extract the edges
242   align = vec_lvsr( 0, to );                                // permute map to misalign data
243   MSQ = vec_perm(edges,(__vector unsigned char)from,align);   // misalign the data (MSQ)
244   LSQ = vec_perm((__vector unsigned char)from,edges,align);   // misalign the data (LSQ)
245   vec_st( LSQ, 15, (unsigned char *)to );                   // Store the LSQ part first
246   vec_st( MSQ, 0, (unsigned char *)to );                    // Store the MSQ part
247 }
248
249 template<> inline float  ei_pfirst(const v4f&  a)
250 {
251   float __attribute__(aligned(16)) af[4];
252   vec_st(a, 0, af);
253   return af[0];
254 }
255
256 template<> inline int    ei_pfirst(const v4i&  a)
257 {
258   int __attribute__(aligned(16)) ai[4];
259   vec_st(a, 0, ai);
260   return ai[0];
261 }
262
263 inline v4f ei_preduxp(const v4f* vecs)
264 {
265   v4f v[4], sum[4];
266
267   // It's easier and faster to transpose then add as columns
268   // Check: http://www.freevec.org/function/matrix_4x4_transpose_floats for explanation
269   // Do the transpose, first set of moves
270   v[0] = vec_mergeh(vecs[0], vecs[2]);
271   v[1] = vec_mergel(vecs[0], vecs[2]);
272   v[2] = vec_mergeh(vecs[1], vecs[3]);
273   v[3] = vec_mergel(vecs[1], vecs[3]);
274   // Get the resulting vectors
275   sum[0] = vec_mergeh(v[0], v[2]);
276   sum[1] = vec_mergel(v[0], v[2]);
277   sum[2] = vec_mergeh(v[1], v[3]);
278   sum[3] = vec_mergel(v[1], v[3]);
279
280   // Now do the summation:
281   // Lines 0+1
282   sum[0] = vec_add(sum[0], sum[1]);
283   // Lines 2+3
284   sum[1] = vec_add(sum[2], sum[3]);
285   // Add the results
286   sum[0] = vec_add(sum[0], sum[1]);
287   return sum[0];
288 }
289
290 inline float ei_predux(const v4f& a)
291 {
292   v4f b, sum;
293   b = (v4f)vec_sld(a, a, 8);
294   sum = vec_add(a, b);
295   b = (v4f)vec_sld(sum, sum, 4);
296   sum = vec_add(sum, b);
297   return ei_pfirst(sum);
298 }
299
300 inline v4i  ei_preduxp(const v4i* vecs)
301 {
302   v4i v[4], sum[4];
303
304   // It's easier and faster to transpose then add as columns
305   // Check: http://www.freevec.org/function/matrix_4x4_transpose_floats for explanation
306   // Do the transpose, first set of moves
307   v[0] = vec_mergeh(vecs[0], vecs[2]);
308   v[1] = vec_mergel(vecs[0], vecs[2]);
309   v[2] = vec_mergeh(vecs[1], vecs[3]);
310   v[3] = vec_mergel(vecs[1], vecs[3]);
311   // Get the resulting vectors
312   sum[0] = vec_mergeh(v[0], v[2]);
313   sum[1] = vec_mergel(v[0], v[2]);
314   sum[2] = vec_mergeh(v[1], v[3]);
315   sum[3] = vec_mergel(v[1], v[3]);
316
317   // Now do the summation:
318   // Lines 0+1
319   sum[0] = vec_add(sum[0], sum[1]);
320   // Lines 2+3
321   sum[1] = vec_add(sum[2], sum[3]);
322   // Add the results
323   sum[0] = vec_add(sum[0], sum[1]);
324   return sum[0];
325 }
326
327 inline int ei_predux(const v4i& a)
328 {
329   USE_CONST_v0i;
330   v4i sum;
331   sum = vec_sums(a, v0i);
332   sum = vec_sld(sum, v0i, 12);
333   return ei_pfirst(sum);
334 }
335
336 template<int Offset>
337 struct ei_palign_impl<Offset, v4f>
338 {
339   inline static void run(v4f& first, const v4f& second)
340   {
341     first = vec_sld(first, second, Offset*4);
342   }
343 };
344
345 template<int Offset>
346 struct ei_palign_impl<Offset, v4i>
347 {
348   inline static void run(v4i& first, const v4i& second)
349   {
350     first = vec_sld(first, second, Offset*4);
351   }
352 };
353
354 #endif // EIGEN_PACKET_MATH_ALTIVEC_H