Mercurial > hg > octave-lyh

comparison liboctave/mx-inlines.cc @ 10481:e8811e5dd699

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
avoid exception throwing in mx-inline loops
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 01 Apr 2010 08:27:23 +0200
parents 532802559f39
children 2645a6b1027b
comparison
equal deleted inserted replaced
10480:19e1e4470e01 10481:e8811e5dd699
39 #include "Array-util.h" 39 #include "Array-util.h"
40 40
41 // Provides some commonly repeated, basic loop templates. 41 // Provides some commonly repeated, basic loop templates.
42 42
43 template <class R, class S> 43 template <class R, class S>
44 inline void mx_inline_fill (size_t n, R *r, S s) 44 inline void mx_inline_fill (size_t n, R *r, S s) throw ()
45 { for (size_t i = 0; i < n; i++) r[i] = s; } 45 { for (size_t i = 0; i < n; i++) r[i] = s; }
46 46
47 #define DEFMXUNOP(F, OP) \ 47 #define DEFMXUNOP(F, OP) \
48 template <class R, class X> \ 48 template <class R, class X> \
49 inline void F (size_t n, R *r, const X *x) \ 49 inline void F (size_t n, R *r, const X *x) throw () \
50 { for (size_t i = 0; i < n; i++) r[i] = OP x[i]; } 50 { for (size_t i = 0; i < n; i++) r[i] = OP x[i]; }
51 51
52 DEFMXUNOP (mx_inline_uminus, -) 52 DEFMXUNOP (mx_inline_uminus, -)
53 53
54 #define DEFMXUNOPEQ(F, OP) \ 54 #define DEFMXUNOPEQ(F, OP) \
55 template <class R> \ 55 template <class R> \
56 inline void F (size_t n, R *r) \ 56 inline void F (size_t n, R *r) throw () \
57 { for (size_t i = 0; i < n; i++) r[i] = OP r[i]; } 57 { for (size_t i = 0; i < n; i++) r[i] = OP r[i]; }
58 58
59 DEFMXUNOPEQ (mx_inline_uminus2, -) 59 DEFMXUNOPEQ (mx_inline_uminus2, -)
60 60
61 #define DEFMXUNBOOLOP(F, OP) \ 61 #define DEFMXUNBOOLOP(F, OP) \
62 template <class X> \ 62 template <class X> \
63 inline void F (size_t n, bool *r, const X *x) \ 63 inline void F (size_t n, bool *r, const X *x) throw () \
64 { const X zero = X(); for (size_t i = 0; i < n; i++) r[i] = x[i] OP zero; } 64 { const X zero = X(); for (size_t i = 0; i < n; i++) r[i] = x[i] OP zero; }
65 65
66 DEFMXUNBOOLOP (mx_inline_iszero, ==) 66 DEFMXUNBOOLOP (mx_inline_iszero, ==)
67 DEFMXUNBOOLOP (mx_inline_notzero, !=) 67 DEFMXUNBOOLOP (mx_inline_notzero, !=)
68 68
69 #define DEFMXBINOP(F, OP) \ 69 #define DEFMXBINOP(F, OP) \
70 template <class R, class X, class Y> \ 70 template <class R, class X, class Y> \
71 inline void F (size_t n, R *r, const X *x, const Y *y) \ 71 inline void F (size_t n, R *r, const X *x, const Y *y) throw () \
72 { for (size_t i = 0; i < n; i++) r[i] = x[i] OP y[i]; } \ 72 { for (size_t i = 0; i < n; i++) r[i] = x[i] OP y[i]; } \
73 template <class R, class X, class Y> \ 73 template <class R, class X, class Y> \
74 inline void F (size_t n, R *r, const X *x, Y y) \ 74 inline void F (size_t n, R *r, const X *x, Y y) throw () \
75 { for (size_t i = 0; i < n; i++) r[i] = x[i] OP y; } \ 75 { for (size_t i = 0; i < n; i++) r[i] = x[i] OP y; } \
76 template <class R, class X, class Y> \ 76 template <class R, class X, class Y> \
77 inline void F (size_t n, R *r, X x, const Y *y) \ 77 inline void F (size_t n, R *r, X x, const Y *y) throw () \
78 { for (size_t i = 0; i < n; i++) r[i] = x OP y[i]; } 78 { for (size_t i = 0; i < n; i++) r[i] = x OP y[i]; }
79 79
80 DEFMXBINOP (mx_inline_add, +) 80 DEFMXBINOP (mx_inline_add, +)
81 DEFMXBINOP (mx_inline_sub, -) 81 DEFMXBINOP (mx_inline_sub, -)
82 DEFMXBINOP (mx_inline_mul, *) 82 DEFMXBINOP (mx_inline_mul, *)
83 DEFMXBINOP (mx_inline_div, /) 83 DEFMXBINOP (mx_inline_div, /)
84 84
85 #define DEFMXBINOPEQ(F, OP) \ 85 #define DEFMXBINOPEQ(F, OP) \
86 template <class R, class X> \ 86 template <class R, class X> \
87 inline void F (size_t n, R *r, const X *x) \ 87 inline void F (size_t n, R *r, const X *x) throw () \
88 { for (size_t i = 0; i < n; i++) r[i] OP x[i]; } \ 88 { for (size_t i = 0; i < n; i++) r[i] OP x[i]; } \
89 template <class R, class X> \ 89 template <class R, class X> \
90 inline void F (size_t n, R *r, X x) \ 90 inline void F (size_t n, R *r, X x) throw () \
91 { for (size_t i = 0; i < n; i++) r[i] OP x; } 91 { for (size_t i = 0; i < n; i++) r[i] OP x; }
92 92
93 DEFMXBINOPEQ (mx_inline_add2, +=) 93 DEFMXBINOPEQ (mx_inline_add2, +=)
94 DEFMXBINOPEQ (mx_inline_sub2, -=) 94 DEFMXBINOPEQ (mx_inline_sub2, -=)
95 DEFMXBINOPEQ (mx_inline_mul2, *=) 95 DEFMXBINOPEQ (mx_inline_mul2, *=)
96 DEFMXBINOPEQ (mx_inline_div2, /=) 96 DEFMXBINOPEQ (mx_inline_div2, /=)
97 97
98 #define DEFMXCMPOP(F, OP) \ 98 #define DEFMXCMPOP(F, OP) \
99 template <class X, class Y> \ 99 template <class X, class Y> \
100 inline void F (size_t n, bool *r, const X *x, const Y *y) \ 100 inline void F (size_t n, bool *r, const X *x, const Y *y) throw () \
101 { for (size_t i = 0; i < n; i++) r[i] = x[i] OP y[i]; } \ 101 { for (size_t i = 0; i < n; i++) r[i] = x[i] OP y[i]; } \
102 template <class X, class Y> \ 102 template <class X, class Y> \
103 inline void F (size_t n, bool *r, const X *x, Y y) \ 103 inline void F (size_t n, bool *r, const X *x, Y y) throw () \
104 { for (size_t i = 0; i < n; i++) r[i] = x[i] OP y; } \ 104 { for (size_t i = 0; i < n; i++) r[i] = x[i] OP y; } \
105 template <class X, class Y> \ 105 template <class X, class Y> \
106 inline void F (size_t n, bool *r, X x, const Y *y) \ 106 inline void F (size_t n, bool *r, X x, const Y *y) throw () \
107 { for (size_t i = 0; i < n; i++) r[i] = x OP y[i]; } 107 { for (size_t i = 0; i < n; i++) r[i] = x OP y[i]; }
108 108
109 DEFMXCMPOP (mx_inline_lt, <) 109 DEFMXCMPOP (mx_inline_lt, <)
110 DEFMXCMPOP (mx_inline_le, <=) 110 DEFMXCMPOP (mx_inline_le, <=)
111 DEFMXCMPOP (mx_inline_gt, >) 111 DEFMXCMPOP (mx_inline_gt, >)
113 DEFMXCMPOP (mx_inline_eq, ==) 113 DEFMXCMPOP (mx_inline_eq, ==)
114 DEFMXCMPOP (mx_inline_ne, !=) 114 DEFMXCMPOP (mx_inline_ne, !=)
115 115
116 // Convert to logical value, for logical op purposes. 116 // Convert to logical value, for logical op purposes.
117 template <class T> inline bool logical_value (T x) { return x; } 117 template <class T> inline bool logical_value (T x) { return x; }
118 template <class T> inline bool logical_value (const std::complex<T>& x)
119 { return x.real () != 0 && x.imag () != 0; }
118 template <class T> inline bool logical_value (const octave_int<T>& x) 120 template <class T> inline bool logical_value (const octave_int<T>& x)
119 { return x.value (); } 121 { return x.value (); }
120 122
121 // NaNs in real data should generate an error. Doing it on-the-fly is faster.
122
123 #define DEFLOGCHKNAN(ARG, ZERO) \
124 inline bool logical_value (ARG x) \
125 { if (xisnan (x)) gripe_nan_to_logical_conversion (); return x != ZERO; }
126
127 DEFLOGCHKNAN (double, 0.0)
128 DEFLOGCHKNAN (const Complex&, 0.0)
129 DEFLOGCHKNAN (float, 0.0f)
130 DEFLOGCHKNAN (const FloatComplex&, 0.0f)
131
132 template <class X> 123 template <class X>
133 void mx_inline_not (size_t n, bool *r, const X* x) 124 void mx_inline_not (size_t n, bool *r, const X* x) throw ()
134 { 125 {
135 for (size_t i = 0; i < n; i++) 126 for (size_t i = 0; i < n; i++)
136 r[i] = ! logical_value (x[i]); 127 r[i] = ! logical_value (x[i]);
137 } 128 }
138 129
139 inline void mx_inline_not2 (size_t n, bool *r) 130 inline void mx_inline_not2 (size_t n, bool *r) throw ()
140 { 131 {
141 for (size_t i = 0; i < n; i++) r[i] = ! r[i]; 132 for (size_t i = 0; i < n; i++) r[i] = ! r[i];
142 } 133 }
143 134
144 #define DEFMXBOOLOP(F, NOT1, OP, NOT2) \ 135 #define DEFMXBOOLOP(F, NOT1, OP, NOT2) \
145 template <class X, class Y> \ 136 template <class X, class Y> \
146 inline void F (size_t n, bool *r, const X *x, const Y *y) \ 137 inline void F (size_t n, bool *r, const X *x, const Y *y) throw () \
147 { \ 138 { \
148 for (size_t i = 0; i < n; i++) \ 139 for (size_t i = 0; i < n; i++) \
149 r[i] = (NOT1 logical_value (x[i])) OP (NOT2 logical_value (y[i])); \ 140 r[i] = (NOT1 logical_value (x[i])) OP (NOT2 logical_value (y[i])); \
150 } \ 141 } \
151 template <class X, class Y> \ 142 template <class X, class Y> \
152 inline void F (size_t n, bool *r, const X *x, Y y) \ 143 inline void F (size_t n, bool *r, const X *x, Y y) throw () \
153 { \ 144 { \
154 const bool yy = (NOT2 logical_value (y)); \ 145 const bool yy = (NOT2 logical_value (y)); \
155 for (size_t i = 0; i < n; i++) \ 146 for (size_t i = 0; i < n; i++) \
156 r[i] = (NOT1 logical_value (x[i])) OP yy; \ 147 r[i] = (NOT1 logical_value (x[i])) OP yy; \
157 } \ 148 } \
158 template <class X, class Y> \ 149 template <class X, class Y> \
159 inline void F (size_t n, bool *r, X x, const Y *y) \ 150 inline void F (size_t n, bool *r, X x, const Y *y) throw () \
160 { \ 151 { \
161 const bool xx = (NOT1 logical_value (x)); \ 152 const bool xx = (NOT1 logical_value (x)); \
162 for (size_t i = 0; i < n; i++) \ 153 for (size_t i = 0; i < n; i++) \
163 r[i] = xx OP (NOT2 logical_value (y[i])); \ 154 r[i] = xx OP (NOT2 logical_value (y[i])); \
164 } 155 }
170 DEFMXBOOLOP (mx_inline_and_not, , &, !) 161 DEFMXBOOLOP (mx_inline_and_not, , &, !)
171 DEFMXBOOLOP (mx_inline_or_not, , |, !) 162 DEFMXBOOLOP (mx_inline_or_not, , |, !)
172 163
173 #define DEFMXBOOLOPEQ(F, OP) \ 164 #define DEFMXBOOLOPEQ(F, OP) \
174 template <class X> \ 165 template <class X> \
175 inline void F (size_t n, bool *r, const X *x) \ 166 inline void F (size_t n, bool *r, const X *x) throw () \
176 { \ 167 { \
177 for (size_t i = 0; i < n; i++) \ 168 for (size_t i = 0; i < n; i++) \
178 r[i] OP logical_value (x[i]); \ 169 r[i] OP logical_value (x[i]); \
179 } \ 170 } \
180 171
181 DEFMXBOOLOPEQ (mx_inline_and2, &=) 172 DEFMXBOOLOPEQ (mx_inline_and2, &=)
182 DEFMXBOOLOPEQ (mx_inline_or2, |=) 173 DEFMXBOOLOPEQ (mx_inline_or2, |=)
183 174
184 template <class T> 175 template <class T>
185 inline bool 176 inline bool
186 mx_inline_any_nan (size_t n, const T* x) 177 mx_inline_any_nan (size_t n, const T* x) throw ()
187 { 178 {
188 for (size_t i = 0; i < n; i++) 179 for (size_t i = 0; i < n; i++)
189 { 180 {
190 if (xisnan (x[i])) 181 if (xisnan (x[i]))
191 return true; 182 return true;
194 return false; 185 return false;
195 } 186 }
196 187
197 template <class T> 188 template <class T>
198 inline bool 189 inline bool
199 mx_inline_any_negative (size_t n, const T* x) 190 mx_inline_any_negative (size_t n, const T* x) throw ()
200 { 191 {
201 for (size_t i = 0; i < n; i++) 192 for (size_t i = 0; i < n; i++)
202 { 193 {
203 if (x[i] < 0) 194 if (x[i] < 0)
204 return true; 195 return true;
207 return false; 198 return false;
208 } 199 }
209 200
210 template<class T> 201 template<class T>
211 inline bool 202 inline bool
212 mx_inline_all_real (size_t n, const std::complex<T>* x) 203 mx_inline_all_real (size_t n, const std::complex<T>* x) throw ()
213 { 204 {
214 for (size_t i = 0; i < n; i++) 205 for (size_t i = 0; i < n; i++)
215 { 206 {
216 if (x[i].imag () != 0) 207 if (x[i].imag () != 0)
217 return false; 208 return false;
220 return true; 211 return true;
221 } 212 }
222 213
223 #define DEFMXMAPPER(F, FUN) \ 214 #define DEFMXMAPPER(F, FUN) \
224 template <class T> \ 215 template <class T> \
225 inline void F (size_t n, T *r, const T *x) \ 216 inline void F (size_t n, T *r, const T *x) throw () \
226 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i]); } 217 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i]); }
227 218
228 template<class T> 219 template<class T>
229 inline void mx_inline_real (size_t n, T *r, const std::complex<T>* x) 220 inline void mx_inline_real (size_t n, T *r, const std::complex<T>* x) throw ()
230 { for (size_t i = 0; i < n; i++) r[i] = x[i].real (); } 221 { for (size_t i = 0; i < n; i++) r[i] = x[i].real (); }
231 template<class T> 222 template<class T>
232 inline void mx_inline_imag (size_t n, T *r, const std::complex<T>* x) 223 inline void mx_inline_imag (size_t n, T *r, const std::complex<T>* x) throw ()
233 { for (size_t i = 0; i < n; i++) r[i] = x[i].imag (); } 224 { for (size_t i = 0; i < n; i++) r[i] = x[i].imag (); }
234 225
235 // Pairwise minimums/maximums 226 // Pairwise minimums/maximums
236 #define DEFMXMAPPER2(F, FUN) \ 227 #define DEFMXMAPPER2(F, FUN) \
237 template <class T> \ 228 template <class T> \
238 inline void F (size_t n, T *r, const T *x, const T *y) \ 229 inline void F (size_t n, T *r, const T *x, const T *y) throw () \
239 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i], y[i]); } \ 230 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i], y[i]); } \
240 template <class T> \ 231 template <class T> \
241 inline void F (size_t n, T *r, const T *x, T y) \ 232 inline void F (size_t n, T *r, const T *x, T y) throw () \
242 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i], y); } \ 233 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i], y); } \
243 template <class T> \ 234 template <class T> \
244 inline void F (size_t n, T *r, T x, const T *y) \ 235 inline void F (size_t n, T *r, T x, const T *y) throw () \
245 { for (size_t i = 0; i < n; i++) r[i] = FUN (x, y[i]); } 236 { for (size_t i = 0; i < n; i++) r[i] = FUN (x, y[i]); }
246 237
247 DEFMXMAPPER2 (mx_inline_xmin, xmin) 238 DEFMXMAPPER2 (mx_inline_xmin, xmin)
248 DEFMXMAPPER2 (mx_inline_xmax, xmax) 239 DEFMXMAPPER2 (mx_inline_xmax, xmax)
249 240
250 // Specialize array-scalar max/min 241 // Specialize array-scalar max/min
251 #define DEFMINMAXSPEC(T, F, OP) \ 242 #define DEFMINMAXSPEC(T, F, OP) \
252 template <> \ 243 template <> \
253 inline void F<T> (size_t n, T *r, const T *x, T y) \ 244 inline void F<T> (size_t n, T *r, const T *x, T y) throw () \
254 { \ 245 { \
255 if (xisnan (y)) \ 246 if (xisnan (y)) \
256 std::memcpy (r, x, n * sizeof (T)); \ 247 std::memcpy (r, x, n * sizeof (T)); \
257 else \ 248 else \
258 for (size_t i = 0; i < n; i++) r[i] = (x[i] OP y) ? x[i] : y; \ 249 for (size_t i = 0; i < n; i++) r[i] = (x[i] OP y) ? x[i] : y; \
259 } \ 250 } \
260 template <> \ 251 template <> \
261 inline void F<T> (size_t n, T *r, T x, const T *y) \ 252 inline void F<T> (size_t n, T *r, T x, const T *y) throw () \
262 { \ 253 { \
263 if (xisnan (x)) \ 254 if (xisnan (x)) \
264 std::memcpy (r, y, n * sizeof (T)); \ 255 std::memcpy (r, y, n * sizeof (T)); \
265 else \ 256 else \
266 for (size_t i = 0; i < n; i++) r[i] = (y[i] OP x) ? y[i] : x; \ 257 for (size_t i = 0; i < n; i++) r[i] = (y[i] OP x) ? y[i] : x; \
272 DEFMINMAXSPEC (float, mx_inline_xmax, >=) 263 DEFMINMAXSPEC (float, mx_inline_xmax, >=)
273 264
274 // Pairwise power 265 // Pairwise power
275 #define DEFMXMAPPER2X(F, FUN) \ 266 #define DEFMXMAPPER2X(F, FUN) \
276 template <class R, class X, class Y> \ 267 template <class R, class X, class Y> \
277 inline void F (size_t n, R *r, const X *x, const Y *y) \ 268 inline void F (size_t n, R *r, const X *x, const Y *y) throw () \
278 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i], y[i]); } \ 269 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i], y[i]); } \
279 template <class R, class X, class Y> \ 270 template <class R, class X, class Y> \
280 inline void F (size_t n, R *r, const X *x, Y y) \ 271 inline void F (size_t n, R *r, const X *x, Y y) throw () \
281 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i], y); } \ 272 { for (size_t i = 0; i < n; i++) r[i] = FUN (x[i], y); } \
282 template <class R, class X, class Y> \ 273 template <class R, class X, class Y> \
283 inline void F (size_t n, R *r, X x, const Y *y) \ 274 inline void F (size_t n, R *r, X x, const Y *y) throw () \
284 { for (size_t i = 0; i < n; i++) r[i] = FUN (x, y[i]); } 275 { for (size_t i = 0; i < n; i++) r[i] = FUN (x, y[i]); }
285 276
286 DEFMXMAPPER2X (mx_inline_pow, std::pow) 277 DEFMXMAPPER2X (mx_inline_pow, std::pow)
287 278
288 // Arbitrary function appliers. The function is a template parameter to enable 279 // Arbitrary function appliers. The function is a template parameter to enable
289 // inlining. 280 // inlining.
290 template <class R, class X, R fun (X x)> 281 template <class R, class X, R fun (X x)>
291 inline void mx_inline_map (size_t n, R *r, const X *x) 282 inline void mx_inline_map (size_t n, R *r, const X *x) throw ()
292 { for (size_t i = 0; i < n; i++) r[i] = fun (x[i]); } 283 { for (size_t i = 0; i < n; i++) r[i] = fun (x[i]); }
293 284
294 template <class R, class X, R fun (const X& x)> 285 template <class R, class X, R fun (const X& x)>
295 inline void mx_inline_map (size_t n, R *r, const X *x) 286 inline void mx_inline_map (size_t n, R *r, const X *x) throw ()
296 { for (size_t i = 0; i < n; i++) r[i] = fun (x[i]); } 287 { for (size_t i = 0; i < n; i++) r[i] = fun (x[i]); }
297 288
298 // Appliers. Since these call the operation just once, we pass it as 289 // Appliers. Since these call the operation just once, we pass it as
299 // a pointer, to allow the compiler reduce number of instances. 290 // a pointer, to allow the compiler reduce number of instances.
300 291
301 template <class R, class X> 292 template <class R, class X>
302 inline Array<R> 293 inline Array<R>
303 do_mx_unary_op (const Array<X>& x, 294 do_mx_unary_op (const Array<X>& x,
304 void (*op) (size_t, R *, const X *)) 295 void (*op) (size_t, R *, const X *) throw ())
305 { 296 {
306 Array<R> r (x.dims ()); 297 Array<R> r (x.dims ());
307 op (r.numel (), r.fortran_vec (), x.data ()); 298 op (r.numel (), r.fortran_vec (), x.data ());
308 return r; 299 return r;
309 } 300 }
325 } 316 }
326 317
327 template <class R> 318 template <class R>
328 inline Array<R>& 319 inline Array<R>&
329 do_mx_inplace_op (Array<R>& r, 320 do_mx_inplace_op (Array<R>& r,
330 void (*op) (size_t, R *)) 321 void (*op) (size_t, R *) throw ())
331 { 322 {
332 op (r.numel (), r.fortran_vec ()); 323 op (r.numel (), r.fortran_vec ());
333 return r; 324 return r;
334 } 325 }
335 326
336 327
337 template <class R, class X, class Y> 328 template <class R, class X, class Y>
338 inline Array<R> 329 inline Array<R>
339 do_mm_binary_op (const Array<X>& x, const Array<Y>& y, 330 do_mm_binary_op (const Array<X>& x, const Array<Y>& y,
340 void (*op) (size_t, R *, const X *, const Y *), 331 void (*op) (size_t, R *, const X *, const Y *) throw (),
341 const char *opname) 332 const char *opname)
342 { 333 {
343 dim_vector dx = x.dims (), dy = y.dims (); 334 dim_vector dx = x.dims (), dy = y.dims ();
344 if (dx == dy) 335 if (dx == dy)
345 { 336 {
355 } 346 }
356 347
357 template <class R, class X, class Y> 348 template <class R, class X, class Y>
358 inline Array<R> 349 inline Array<R>
359 do_ms_binary_op (const Array<X>& x, const Y& y, 350 do_ms_binary_op (const Array<X>& x, const Y& y,
360 void (*op) (size_t, R *, const X *, Y)) 351 void (*op) (size_t, R *, const X *, Y) throw ())
361 { 352 {
362 Array<R> r (x.dims ()); 353 Array<R> r (x.dims ());
363 op (r.length (), r.fortran_vec (), x.data (), y); 354 op (r.length (), r.fortran_vec (), x.data (), y);
364 return r; 355 return r;
365 } 356 }
366 357
367 template <class R, class X, class Y> 358 template <class R, class X, class Y>
368 inline Array<R> 359 inline Array<R>
369 do_sm_binary_op (const X& x, const Array<Y>& y, 360 do_sm_binary_op (const X& x, const Array<Y>& y,
370 void (*op) (size_t, R *, X, const Y *)) 361 void (*op) (size_t, R *, X, const Y *) throw ())
371 { 362 {
372 Array<R> r (y.dims ()); 363 Array<R> r (y.dims ());
373 op (r.length (), r.fortran_vec (), x, y.data ()); 364 op (r.length (), r.fortran_vec (), x, y.data ());
374 return r; 365 return r;
375 } 366 }
376 367
377 template <class R, class X> 368 template <class R, class X>
378 inline Array<R>& 369 inline Array<R>&
379 do_mm_inplace_op (Array<R>& r, const Array<X>& x, 370 do_mm_inplace_op (Array<R>& r, const Array<X>& x,
380 void (*op) (size_t, R *, const X *), 371 void (*op) (size_t, R *, const X *) throw (),
381 const char *opname) 372 const char *opname)
382 { 373 {
383 dim_vector dr = r.dims (), dx = x.dims (); 374 dim_vector dr = r.dims (), dx = x.dims ();
384 if (dr == dx) 375 if (dr == dx)
385 op (r.length (), r.fortran_vec (), x.data ()); 376 op (r.length (), r.fortran_vec (), x.data ());
389 } 380 }
390 381
391 template <class R, class X> 382 template <class R, class X>
392 inline Array<R>& 383 inline Array<R>&
393 do_ms_inplace_op (Array<R>& r, const X& x, 384 do_ms_inplace_op (Array<R>& r, const X& x,
394 void (*op) (size_t, R *, X)) 385 void (*op) (size_t, R *, X) throw ())
395 { 386 {
396 op (r.length (), r.fortran_vec (), x); 387 op (r.length (), r.fortran_vec (), x);
397 return r; 388 return r;
398 } 389 }
399 390
400 template <class T1, class T2> 391 template <class T1, class T2>
401 inline bool 392 inline bool
402 mx_inline_equal (size_t n, const T1 *x, const T2 *y) 393 mx_inline_equal (size_t n, const T1 *x, const T2 *y) throw ()
403 { 394 {
404 for (size_t i = 0; i < n; i++) 395 for (size_t i = 0; i < n; i++)
405 if (x[i] != y[i]) 396 if (x[i] != y[i])
406 return false; 397 return false;
407 return true; 398 return true;
408 } 399 }
409 400
410 // FIXME: Due to a performance defect in g++ (<= 4.3), std::norm is slow unless 401 template <class T>
411 // ffast-math is on (not by default even with -O3). The following helper function 402 inline bool
412 // gives the expected straightforward implementation of std::norm. 403 do_mx_check (const Array<T>& a,
404 bool (*op) (size_t, const T *) throw ())
405 {
406 return op (a.numel (), a.data ());
407 }
408
409 // NOTE: we don't use std::norm because it typically does some heavyweight
410 // magic to avoid underflows, which we don't need here.
413 template <class T> 411 template <class T>
414 inline T cabsq (const std::complex<T>& c) 412 inline T cabsq (const std::complex<T>& c)
415 { return c.real () * c.real () + c.imag () * c.imag (); } 413 { return c.real () * c.real () + c.imag () * c.imag (); }
416 414
417 // default. works for integers and bool. 415 // default. works for integers and bool.