]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // |
2 | // Copyright (c) 2000-2009 | |
3 | // Joerg Walter, Mathias Koch, Gunter Winkler | |
4 | // | |
5 | // Distributed under the Boost Software License, Version 1.0. (See | |
6 | // accompanying file LICENSE_1_0.txt or copy at | |
7 | // http://www.boost.org/LICENSE_1_0.txt) | |
8 | // | |
9 | // The authors gratefully acknowledge the support of | |
10 | // GeNeSys mbH & Co. KG in producing this work. | |
11 | // | |
12 | ||
13 | #ifndef _BOOST_UBLAS_FUNCTIONAL_ | |
14 | #define _BOOST_UBLAS_FUNCTIONAL_ | |
15 | ||
16 | #include <functional> | |
17 | ||
18 | #include <boost/core/ignore_unused.hpp> | |
19 | ||
20 | #include <boost/numeric/ublas/traits.hpp> | |
21 | #ifdef BOOST_UBLAS_USE_DUFF_DEVICE | |
22 | #include <boost/numeric/ublas/detail/duff.hpp> | |
23 | #endif | |
24 | #ifdef BOOST_UBLAS_USE_SIMD | |
25 | #include <boost/numeric/ublas/detail/raw.hpp> | |
26 | #else | |
27 | namespace boost { namespace numeric { namespace ublas { namespace raw { | |
28 | }}}} | |
29 | #endif | |
30 | #ifdef BOOST_UBLAS_HAVE_BINDINGS | |
31 | #include <boost/numeric/bindings/traits/std_vector.hpp> | |
32 | #include <boost/numeric/bindings/traits/ublas_vector.hpp> | |
33 | #include <boost/numeric/bindings/traits/ublas_matrix.hpp> | |
34 | #include <boost/numeric/bindings/atlas/cblas.hpp> | |
35 | #endif | |
36 | ||
37 | #include <boost/numeric/ublas/detail/definitions.hpp> | |
38 | ||
39 | ||
40 | ||
41 | namespace boost { namespace numeric { namespace ublas { | |
42 | ||
43 | // Scalar functors | |
44 | ||
45 | // Unary | |
46 | template<class T> | |
47 | struct scalar_unary_functor { | |
48 | typedef T value_type; | |
49 | typedef typename type_traits<T>::const_reference argument_type; | |
50 | typedef typename type_traits<T>::value_type result_type; | |
51 | }; | |
52 | ||
53 | template<class T> | |
54 | struct scalar_identity: | |
55 | public scalar_unary_functor<T> { | |
56 | typedef typename scalar_unary_functor<T>::argument_type argument_type; | |
57 | typedef typename scalar_unary_functor<T>::result_type result_type; | |
58 | ||
59 | static BOOST_UBLAS_INLINE | |
60 | result_type apply (argument_type t) { | |
61 | return t; | |
62 | } | |
63 | }; | |
64 | template<class T> | |
65 | struct scalar_negate: | |
66 | public scalar_unary_functor<T> { | |
67 | typedef typename scalar_unary_functor<T>::argument_type argument_type; | |
68 | typedef typename scalar_unary_functor<T>::result_type result_type; | |
69 | ||
70 | static BOOST_UBLAS_INLINE | |
71 | result_type apply (argument_type t) { | |
72 | return - t; | |
73 | } | |
74 | }; | |
75 | template<class T> | |
76 | struct scalar_conj: | |
77 | public scalar_unary_functor<T> { | |
78 | typedef typename scalar_unary_functor<T>::value_type value_type; | |
79 | typedef typename scalar_unary_functor<T>::argument_type argument_type; | |
80 | typedef typename scalar_unary_functor<T>::result_type result_type; | |
81 | ||
82 | static BOOST_UBLAS_INLINE | |
83 | result_type apply (argument_type t) { | |
84 | return type_traits<value_type>::conj (t); | |
85 | } | |
86 | }; | |
87 | ||
88 | // Unary returning real | |
89 | template<class T> | |
90 | struct scalar_real_unary_functor { | |
91 | typedef T value_type; | |
92 | typedef typename type_traits<T>::const_reference argument_type; | |
93 | typedef typename type_traits<T>::real_type result_type; | |
94 | }; | |
95 | ||
96 | template<class T> | |
97 | struct scalar_real: | |
98 | public scalar_real_unary_functor<T> { | |
99 | typedef typename scalar_real_unary_functor<T>::value_type value_type; | |
100 | typedef typename scalar_real_unary_functor<T>::argument_type argument_type; | |
101 | typedef typename scalar_real_unary_functor<T>::result_type result_type; | |
102 | ||
103 | static BOOST_UBLAS_INLINE | |
104 | result_type apply (argument_type t) { | |
105 | return type_traits<value_type>::real (t); | |
106 | } | |
107 | }; | |
108 | template<class T> | |
109 | struct scalar_imag: | |
110 | public scalar_real_unary_functor<T> { | |
111 | typedef typename scalar_real_unary_functor<T>::value_type value_type; | |
112 | typedef typename scalar_real_unary_functor<T>::argument_type argument_type; | |
113 | typedef typename scalar_real_unary_functor<T>::result_type result_type; | |
114 | ||
115 | static BOOST_UBLAS_INLINE | |
116 | result_type apply (argument_type t) { | |
117 | return type_traits<value_type>::imag (t); | |
118 | } | |
119 | }; | |
120 | ||
121 | // Binary | |
122 | template<class T1, class T2> | |
123 | struct scalar_binary_functor { | |
124 | typedef typename type_traits<T1>::const_reference argument1_type; | |
125 | typedef typename type_traits<T2>::const_reference argument2_type; | |
126 | typedef typename promote_traits<T1, T2>::promote_type result_type; | |
127 | }; | |
128 | ||
129 | template<class T1, class T2> | |
130 | struct scalar_plus: | |
131 | public scalar_binary_functor<T1, T2> { | |
132 | typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; | |
133 | typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; | |
134 | typedef typename scalar_binary_functor<T1, T2>::result_type result_type; | |
135 | ||
136 | static BOOST_UBLAS_INLINE | |
137 | result_type apply (argument1_type t1, argument2_type t2) { | |
138 | return t1 + t2; | |
139 | } | |
140 | }; | |
141 | template<class T1, class T2> | |
142 | struct scalar_minus: | |
143 | public scalar_binary_functor<T1, T2> { | |
144 | typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; | |
145 | typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; | |
146 | typedef typename scalar_binary_functor<T1, T2>::result_type result_type; | |
147 | ||
148 | static BOOST_UBLAS_INLINE | |
149 | result_type apply (argument1_type t1, argument2_type t2) { | |
150 | return t1 - t2; | |
151 | } | |
152 | }; | |
153 | template<class T1, class T2> | |
154 | struct scalar_multiplies: | |
155 | public scalar_binary_functor<T1, T2> { | |
156 | typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; | |
157 | typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; | |
158 | typedef typename scalar_binary_functor<T1, T2>::result_type result_type; | |
159 | ||
160 | static BOOST_UBLAS_INLINE | |
161 | result_type apply (argument1_type t1, argument2_type t2) { | |
162 | return t1 * t2; | |
163 | } | |
164 | }; | |
165 | template<class T1, class T2> | |
166 | struct scalar_divides: | |
167 | public scalar_binary_functor<T1, T2> { | |
168 | typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; | |
169 | typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; | |
170 | typedef typename scalar_binary_functor<T1, T2>::result_type result_type; | |
171 | ||
172 | static BOOST_UBLAS_INLINE | |
173 | result_type apply (argument1_type t1, argument2_type t2) { | |
174 | return t1 / t2; | |
175 | } | |
176 | }; | |
177 | ||
178 | template<class T1, class T2> | |
179 | struct scalar_binary_assign_functor { | |
180 | // ISSUE Remove reference to avoid reference to reference problems | |
181 | typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type; | |
182 | typedef typename type_traits<T2>::const_reference argument2_type; | |
183 | }; | |
184 | ||
185 | struct assign_tag {}; | |
186 | struct computed_assign_tag {}; | |
187 | ||
188 | template<class T1, class T2> | |
189 | struct scalar_assign: | |
190 | public scalar_binary_assign_functor<T1, T2> { | |
191 | typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; | |
192 | typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; | |
193 | #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) | |
194 | static const bool computed ; | |
195 | #else | |
196 | static const bool computed = false ; | |
197 | #endif | |
198 | ||
199 | static BOOST_UBLAS_INLINE | |
200 | void apply (argument1_type t1, argument2_type t2) { | |
201 | t1 = t2; | |
202 | } | |
203 | ||
204 | template<class U1, class U2> | |
205 | struct rebind { | |
206 | typedef scalar_assign<U1, U2> other; | |
207 | }; | |
208 | }; | |
209 | ||
210 | #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) | |
211 | template<class T1, class T2> | |
212 | const bool scalar_assign<T1,T2>::computed = false; | |
213 | #endif | |
214 | ||
215 | template<class T1, class T2> | |
216 | struct scalar_plus_assign: | |
217 | public scalar_binary_assign_functor<T1, T2> { | |
218 | typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; | |
219 | typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; | |
220 | #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) | |
221 | static const bool computed ; | |
222 | #else | |
223 | static const bool computed = true ; | |
224 | #endif | |
225 | ||
226 | static BOOST_UBLAS_INLINE | |
227 | void apply (argument1_type t1, argument2_type t2) { | |
228 | t1 += t2; | |
229 | } | |
230 | ||
231 | template<class U1, class U2> | |
232 | struct rebind { | |
233 | typedef scalar_plus_assign<U1, U2> other; | |
234 | }; | |
235 | }; | |
236 | ||
237 | #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) | |
238 | template<class T1, class T2> | |
239 | const bool scalar_plus_assign<T1,T2>::computed = true; | |
240 | #endif | |
241 | ||
242 | template<class T1, class T2> | |
243 | struct scalar_minus_assign: | |
244 | public scalar_binary_assign_functor<T1, T2> { | |
245 | typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; | |
246 | typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; | |
247 | #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) | |
248 | static const bool computed ; | |
249 | #else | |
250 | static const bool computed = true ; | |
251 | #endif | |
252 | ||
253 | static BOOST_UBLAS_INLINE | |
254 | void apply (argument1_type t1, argument2_type t2) { | |
255 | t1 -= t2; | |
256 | } | |
257 | ||
258 | template<class U1, class U2> | |
259 | struct rebind { | |
260 | typedef scalar_minus_assign<U1, U2> other; | |
261 | }; | |
262 | }; | |
263 | ||
264 | #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) | |
265 | template<class T1, class T2> | |
266 | const bool scalar_minus_assign<T1,T2>::computed = true; | |
267 | #endif | |
268 | ||
269 | template<class T1, class T2> | |
270 | struct scalar_multiplies_assign: | |
271 | public scalar_binary_assign_functor<T1, T2> { | |
272 | typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; | |
273 | typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; | |
274 | static const bool computed = true; | |
275 | ||
276 | static BOOST_UBLAS_INLINE | |
277 | void apply (argument1_type t1, argument2_type t2) { | |
278 | t1 *= t2; | |
279 | } | |
280 | ||
281 | template<class U1, class U2> | |
282 | struct rebind { | |
283 | typedef scalar_multiplies_assign<U1, U2> other; | |
284 | }; | |
285 | }; | |
286 | template<class T1, class T2> | |
287 | struct scalar_divides_assign: | |
288 | public scalar_binary_assign_functor<T1, T2> { | |
289 | typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; | |
290 | typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; | |
291 | static const bool computed ; | |
292 | ||
293 | static BOOST_UBLAS_INLINE | |
294 | void apply (argument1_type t1, argument2_type t2) { | |
295 | t1 /= t2; | |
296 | } | |
297 | ||
298 | template<class U1, class U2> | |
299 | struct rebind { | |
300 | typedef scalar_divides_assign<U1, U2> other; | |
301 | }; | |
302 | }; | |
303 | template<class T1, class T2> | |
304 | const bool scalar_divides_assign<T1,T2>::computed = true; | |
305 | ||
306 | template<class T1, class T2> | |
307 | struct scalar_binary_swap_functor { | |
308 | typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type; | |
309 | typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type; | |
310 | }; | |
311 | ||
312 | template<class T1, class T2> | |
313 | struct scalar_swap: | |
314 | public scalar_binary_swap_functor<T1, T2> { | |
315 | typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type; | |
316 | typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type; | |
317 | ||
318 | static BOOST_UBLAS_INLINE | |
319 | void apply (argument1_type t1, argument2_type t2) { | |
320 | std::swap (t1, t2); | |
321 | } | |
322 | ||
323 | template<class U1, class U2> | |
324 | struct rebind { | |
325 | typedef scalar_swap<U1, U2> other; | |
326 | }; | |
327 | }; | |
328 | ||
329 | // Vector functors | |
330 | ||
331 | // Unary returning scalar | |
332 | template<class V> | |
333 | struct vector_scalar_unary_functor { | |
334 | typedef typename V::value_type value_type; | |
335 | typedef typename V::value_type result_type; | |
336 | }; | |
337 | ||
338 | template<class V> | |
339 | struct vector_sum: | |
340 | public vector_scalar_unary_functor<V> { | |
341 | typedef typename vector_scalar_unary_functor<V>::value_type value_type; | |
342 | typedef typename vector_scalar_unary_functor<V>::result_type result_type; | |
343 | ||
344 | template<class E> | |
345 | static BOOST_UBLAS_INLINE | |
346 | result_type apply (const vector_expression<E> &e) { | |
347 | result_type t = result_type (0); | |
348 | typedef typename E::size_type vector_size_type; | |
349 | vector_size_type size (e ().size ()); | |
350 | for (vector_size_type i = 0; i < size; ++ i) | |
351 | t += e () (i); | |
352 | return t; | |
353 | } | |
354 | // Dense case | |
355 | template<class D, class I> | |
356 | static BOOST_UBLAS_INLINE | |
357 | result_type apply (D size, I it) { | |
358 | result_type t = result_type (0); | |
359 | while (-- size >= 0) | |
360 | t += *it, ++ it; | |
361 | return t; | |
362 | } | |
363 | // Sparse case | |
364 | template<class I> | |
365 | static BOOST_UBLAS_INLINE | |
366 | result_type apply (I it, const I &it_end) { | |
367 | result_type t = result_type (0); | |
368 | while (it != it_end) | |
369 | t += *it, ++ it; | |
370 | return t; | |
371 | } | |
372 | }; | |
373 | ||
374 | // Unary returning real scalar | |
375 | template<class V> | |
376 | struct vector_scalar_real_unary_functor { | |
377 | typedef typename V::value_type value_type; | |
378 | typedef typename type_traits<value_type>::real_type real_type; | |
379 | typedef real_type result_type; | |
380 | }; | |
381 | ||
382 | template<class V> | |
383 | struct vector_norm_1: | |
384 | public vector_scalar_real_unary_functor<V> { | |
385 | typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; | |
386 | typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; | |
387 | typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; | |
388 | ||
389 | template<class E> | |
390 | static BOOST_UBLAS_INLINE | |
391 | result_type apply (const vector_expression<E> &e) { | |
392 | real_type t = real_type (); | |
393 | typedef typename E::size_type vector_size_type; | |
394 | vector_size_type size (e ().size ()); | |
395 | for (vector_size_type i = 0; i < size; ++ i) { | |
396 | real_type u (type_traits<value_type>::type_abs (e () (i))); | |
397 | t += u; | |
398 | } | |
399 | return t; | |
400 | } | |
401 | // Dense case | |
402 | template<class D, class I> | |
403 | static BOOST_UBLAS_INLINE | |
404 | result_type apply (D size, I it) { | |
405 | real_type t = real_type (); | |
406 | while (-- size >= 0) { | |
407 | real_type u (type_traits<value_type>::norm_1 (*it)); | |
408 | t += u; | |
409 | ++ it; | |
410 | } | |
411 | return t; | |
412 | } | |
413 | // Sparse case | |
414 | template<class I> | |
415 | static BOOST_UBLAS_INLINE | |
416 | result_type apply (I it, const I &it_end) { | |
417 | real_type t = real_type (); | |
418 | while (it != it_end) { | |
419 | real_type u (type_traits<value_type>::norm_1 (*it)); | |
420 | t += u; | |
421 | ++ it; | |
422 | } | |
423 | return t; | |
424 | } | |
425 | }; | |
426 | template<class V> | |
427 | struct vector_norm_2: | |
428 | public vector_scalar_real_unary_functor<V> { | |
429 | typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; | |
430 | typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; | |
431 | typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; | |
432 | ||
433 | template<class E> | |
434 | static BOOST_UBLAS_INLINE | |
435 | result_type apply (const vector_expression<E> &e) { | |
436 | #ifndef BOOST_UBLAS_SCALED_NORM | |
437 | real_type t = real_type (); | |
438 | typedef typename E::size_type vector_size_type; | |
439 | vector_size_type size (e ().size ()); | |
440 | for (vector_size_type i = 0; i < size; ++ i) { | |
441 | real_type u (type_traits<value_type>::norm_2 (e () (i))); | |
442 | t += u * u; | |
443 | } | |
444 | return type_traits<real_type>::type_sqrt (t); | |
445 | #else | |
446 | real_type scale = real_type (); | |
447 | real_type sum_squares (1); | |
448 | size_type size (e ().size ()); | |
449 | for (size_type i = 0; i < size; ++ i) { | |
450 | real_type u (type_traits<value_type>::norm_2 (e () (i))); | |
451 | if ( real_type () /* zero */ == u ) continue; | |
452 | if (scale < u) { | |
453 | real_type v (scale / u); | |
454 | sum_squares = sum_squares * v * v + real_type (1); | |
455 | scale = u; | |
456 | } else { | |
457 | real_type v (u / scale); | |
458 | sum_squares += v * v; | |
459 | } | |
460 | } | |
461 | return scale * type_traits<real_type>::type_sqrt (sum_squares); | |
462 | #endif | |
463 | } | |
464 | // Dense case | |
465 | template<class D, class I> | |
466 | static BOOST_UBLAS_INLINE | |
467 | result_type apply (D size, I it) { | |
468 | #ifndef BOOST_UBLAS_SCALED_NORM | |
469 | real_type t = real_type (); | |
470 | while (-- size >= 0) { | |
471 | real_type u (type_traits<value_type>::norm_2 (*it)); | |
472 | t += u * u; | |
473 | ++ it; | |
474 | } | |
475 | return type_traits<real_type>::type_sqrt (t); | |
476 | #else | |
477 | real_type scale = real_type (); | |
478 | real_type sum_squares (1); | |
479 | while (-- size >= 0) { | |
480 | real_type u (type_traits<value_type>::norm_2 (*it)); | |
481 | if (scale < u) { | |
482 | real_type v (scale / u); | |
483 | sum_squares = sum_squares * v * v + real_type (1); | |
484 | scale = u; | |
485 | } else { | |
486 | real_type v (u / scale); | |
487 | sum_squares += v * v; | |
488 | } | |
489 | ++ it; | |
490 | } | |
491 | return scale * type_traits<real_type>::type_sqrt (sum_squares); | |
492 | #endif | |
493 | } | |
494 | // Sparse case | |
495 | template<class I> | |
496 | static BOOST_UBLAS_INLINE | |
497 | result_type apply (I it, const I &it_end) { | |
498 | #ifndef BOOST_UBLAS_SCALED_NORM | |
499 | real_type t = real_type (); | |
500 | while (it != it_end) { | |
501 | real_type u (type_traits<value_type>::norm_2 (*it)); | |
502 | t += u * u; | |
503 | ++ it; | |
504 | } | |
505 | return type_traits<real_type>::type_sqrt (t); | |
506 | #else | |
507 | real_type scale = real_type (); | |
508 | real_type sum_squares (1); | |
509 | while (it != it_end) { | |
510 | real_type u (type_traits<value_type>::norm_2 (*it)); | |
511 | if (scale < u) { | |
512 | real_type v (scale / u); | |
513 | sum_squares = sum_squares * v * v + real_type (1); | |
514 | scale = u; | |
515 | } else { | |
516 | real_type v (u / scale); | |
517 | sum_squares += v * v; | |
518 | } | |
519 | ++ it; | |
520 | } | |
521 | return scale * type_traits<real_type>::type_sqrt (sum_squares); | |
522 | #endif | |
523 | } | |
524 | }; | |
525 | template<class V> | |
526 | struct vector_norm_inf: | |
527 | public vector_scalar_real_unary_functor<V> { | |
528 | typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; | |
529 | typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; | |
530 | typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; | |
531 | ||
532 | template<class E> | |
533 | static BOOST_UBLAS_INLINE | |
534 | result_type apply (const vector_expression<E> &e) { | |
535 | real_type t = real_type (); | |
536 | typedef typename E::size_type vector_size_type; | |
537 | vector_size_type size (e ().size ()); | |
538 | for (vector_size_type i = 0; i < size; ++ i) { | |
539 | real_type u (type_traits<value_type>::norm_inf (e () (i))); | |
540 | if (u > t) | |
541 | t = u; | |
542 | } | |
543 | return t; | |
544 | } | |
545 | // Dense case | |
546 | template<class D, class I> | |
547 | static BOOST_UBLAS_INLINE | |
548 | result_type apply (D size, I it) { | |
549 | real_type t = real_type (); | |
550 | while (-- size >= 0) { | |
551 | real_type u (type_traits<value_type>::norm_inf (*it)); | |
552 | if (u > t) | |
553 | t = u; | |
554 | ++ it; | |
555 | } | |
556 | return t; | |
557 | } | |
558 | // Sparse case | |
559 | template<class I> | |
560 | static BOOST_UBLAS_INLINE | |
561 | result_type apply (I it, const I &it_end) { | |
562 | real_type t = real_type (); | |
563 | while (it != it_end) { | |
564 | real_type u (type_traits<value_type>::norm_inf (*it)); | |
565 | if (u > t) | |
566 | t = u; | |
567 | ++ it; | |
568 | } | |
569 | return t; | |
570 | } | |
571 | }; | |
572 | ||
573 | // Unary returning index | |
574 | template<class V> | |
575 | struct vector_scalar_index_unary_functor { | |
576 | typedef typename V::value_type value_type; | |
577 | typedef typename type_traits<value_type>::real_type real_type; | |
578 | typedef typename V::size_type result_type; | |
579 | }; | |
580 | ||
581 | template<class V> | |
582 | struct vector_index_norm_inf: | |
583 | public vector_scalar_index_unary_functor<V> { | |
584 | typedef typename vector_scalar_index_unary_functor<V>::value_type value_type; | |
585 | typedef typename vector_scalar_index_unary_functor<V>::real_type real_type; | |
586 | typedef typename vector_scalar_index_unary_functor<V>::result_type result_type; | |
587 | ||
588 | template<class E> | |
589 | static BOOST_UBLAS_INLINE | |
590 | result_type apply (const vector_expression<E> &e) { | |
591 | // ISSUE For CBLAS compatibility return 0 index in empty case | |
592 | result_type i_norm_inf (0); | |
593 | real_type t = real_type (); | |
594 | typedef typename E::size_type vector_size_type; | |
595 | vector_size_type size (e ().size ()); | |
596 | for (vector_size_type i = 0; i < size; ++ i) { | |
597 | real_type u (type_traits<value_type>::norm_inf (e () (i))); | |
598 | if (u > t) { | |
599 | i_norm_inf = i; | |
600 | t = u; | |
601 | } | |
602 | } | |
603 | return i_norm_inf; | |
604 | } | |
605 | // Dense case | |
606 | template<class D, class I> | |
607 | static BOOST_UBLAS_INLINE | |
608 | result_type apply (D size, I it) { | |
609 | // ISSUE For CBLAS compatibility return 0 index in empty case | |
610 | result_type i_norm_inf (0); | |
611 | real_type t = real_type (); | |
612 | while (-- size >= 0) { | |
613 | real_type u (type_traits<value_type>::norm_inf (*it)); | |
614 | if (u > t) { | |
615 | i_norm_inf = it.index (); | |
616 | t = u; | |
617 | } | |
618 | ++ it; | |
619 | } | |
620 | return i_norm_inf; | |
621 | } | |
622 | // Sparse case | |
623 | template<class I> | |
624 | static BOOST_UBLAS_INLINE | |
625 | result_type apply (I it, const I &it_end) { | |
626 | // ISSUE For CBLAS compatibility return 0 index in empty case | |
627 | result_type i_norm_inf (0); | |
628 | real_type t = real_type (); | |
629 | while (it != it_end) { | |
630 | real_type u (type_traits<value_type>::norm_inf (*it)); | |
631 | if (u > t) { | |
632 | i_norm_inf = it.index (); | |
633 | t = u; | |
634 | } | |
635 | ++ it; | |
636 | } | |
637 | return i_norm_inf; | |
638 | } | |
639 | }; | |
640 | ||
641 | // Binary returning scalar | |
642 | template<class V1, class V2, class TV> | |
643 | struct vector_scalar_binary_functor { | |
644 | typedef TV value_type; | |
645 | typedef TV result_type; | |
646 | }; | |
647 | ||
648 | template<class V1, class V2, class TV> | |
649 | struct vector_inner_prod: | |
650 | public vector_scalar_binary_functor<V1, V2, TV> { | |
651 | typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type; | |
652 | typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type; | |
653 | ||
654 | template<class C1, class C2> | |
655 | static BOOST_UBLAS_INLINE | |
656 | result_type apply (const vector_container<C1> &c1, | |
657 | const vector_container<C2> &c2) { | |
658 | #ifdef BOOST_UBLAS_USE_SIMD | |
659 | using namespace raw; | |
660 | typedef typename C1::size_type vector_size_type; | |
661 | vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ())); | |
662 | const typename V1::value_type *data1 = data_const (c1 ()); | |
663 | const typename V1::value_type *data2 = data_const (c2 ()); | |
664 | vector_size_type s1 = stride (c1 ()); | |
665 | vector_size_type s2 = stride (c2 ()); | |
666 | result_type t = result_type (0); | |
667 | if (s1 == 1 && s2 == 1) { | |
668 | for (vector_size_type i = 0; i < size; ++ i) | |
669 | t += data1 [i] * data2 [i]; | |
670 | } else if (s2 == 1) { | |
671 | for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1) | |
672 | t += data1 [i1] * data2 [i]; | |
673 | } else if (s1 == 1) { | |
674 | for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2) | |
675 | t += data1 [i] * data2 [i2]; | |
676 | } else { | |
677 | for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2) | |
678 | t += data1 [i1] * data2 [i2]; | |
679 | } | |
680 | return t; | |
681 | #elif defined(BOOST_UBLAS_HAVE_BINDINGS) | |
682 | return boost::numeric::bindings::atlas::dot (c1 (), c2 ()); | |
683 | #else | |
684 | return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2)); | |
685 | #endif | |
686 | } | |
687 | template<class E1, class E2> | |
688 | static BOOST_UBLAS_INLINE | |
689 | result_type apply (const vector_expression<E1> &e1, | |
690 | const vector_expression<E2> &e2) { | |
691 | typedef typename E1::size_type vector_size_type; | |
692 | vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ())); | |
693 | result_type t = result_type (0); | |
694 | #ifndef BOOST_UBLAS_USE_DUFF_DEVICE | |
695 | for (vector_size_type i = 0; i < size; ++ i) | |
696 | t += e1 () (i) * e2 () (i); | |
697 | #else | |
698 | vector_size_type i (0); | |
699 | DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i)); | |
700 | #endif | |
701 | return t; | |
702 | } | |
703 | // Dense case | |
704 | template<class D, class I1, class I2> | |
705 | static BOOST_UBLAS_INLINE | |
706 | result_type apply (D size, I1 it1, I2 it2) { | |
707 | result_type t = result_type (0); | |
708 | #ifndef BOOST_UBLAS_USE_DUFF_DEVICE | |
709 | while (-- size >= 0) | |
710 | t += *it1 * *it2, ++ it1, ++ it2; | |
711 | #else | |
712 | DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); | |
713 | #endif | |
714 | return t; | |
715 | } | |
716 | // Packed case | |
717 | template<class I1, class I2> | |
718 | static BOOST_UBLAS_INLINE | |
719 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { | |
720 | result_type t = result_type (0); | |
721 | typedef typename I1::difference_type vector_difference_type; | |
722 | vector_difference_type it1_size (it1_end - it1); | |
723 | vector_difference_type it2_size (it2_end - it2); | |
724 | vector_difference_type diff (0); | |
725 | if (it1_size > 0 && it2_size > 0) | |
726 | diff = it2.index () - it1.index (); | |
727 | if (diff != 0) { | |
728 | vector_difference_type size = (std::min) (diff, it1_size); | |
729 | if (size > 0) { | |
730 | it1 += size; | |
731 | it1_size -= size; | |
732 | diff -= size; | |
733 | } | |
734 | size = (std::min) (- diff, it2_size); | |
735 | if (size > 0) { | |
736 | it2 += size; | |
737 | it2_size -= size; | |
738 | diff += size; | |
739 | } | |
740 | } | |
741 | vector_difference_type size ((std::min) (it1_size, it2_size)); | |
742 | while (-- size >= 0) | |
743 | t += *it1 * *it2, ++ it1, ++ it2; | |
744 | return t; | |
745 | } | |
746 | // Sparse case | |
747 | template<class I1, class I2> | |
748 | static BOOST_UBLAS_INLINE | |
749 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { | |
750 | result_type t = result_type (0); | |
751 | if (it1 != it1_end && it2 != it2_end) { | |
752 | while (true) { | |
753 | if (it1.index () == it2.index ()) { | |
754 | t += *it1 * *it2, ++ it1, ++ it2; | |
755 | if (it1 == it1_end || it2 == it2_end) | |
756 | break; | |
757 | } else if (it1.index () < it2.index ()) { | |
758 | increment (it1, it1_end, it2.index () - it1.index ()); | |
759 | if (it1 == it1_end) | |
760 | break; | |
761 | } else if (it1.index () > it2.index ()) { | |
762 | increment (it2, it2_end, it1.index () - it2.index ()); | |
763 | if (it2 == it2_end) | |
764 | break; | |
765 | } | |
766 | } | |
767 | } | |
768 | return t; | |
769 | } | |
770 | }; | |
771 | ||
772 | // Matrix functors | |
773 | ||
774 | // Binary returning vector | |
775 | template<class M1, class M2, class TV> | |
776 | struct matrix_vector_binary_functor { | |
777 | typedef typename M1::size_type size_type; | |
778 | typedef typename M1::difference_type difference_type; | |
779 | typedef TV value_type; | |
780 | typedef TV result_type; | |
781 | }; | |
782 | ||
783 | template<class M1, class M2, class TV> | |
784 | struct matrix_vector_prod1: | |
785 | public matrix_vector_binary_functor<M1, M2, TV> { | |
786 | typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type; | |
787 | typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type; | |
788 | typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type; | |
789 | typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type; | |
790 | ||
791 | template<class C1, class C2> | |
792 | static BOOST_UBLAS_INLINE | |
793 | result_type apply (const matrix_container<C1> &c1, | |
794 | const vector_container<C2> &c2, | |
795 | size_type i) { | |
796 | #ifdef BOOST_UBLAS_USE_SIMD | |
797 | using namespace raw; | |
798 | size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ()); | |
799 | const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ()); | |
800 | const typename M2::value_type *data2 = data_const (c2 ()); | |
801 | size_type s1 = stride2 (c1 ()); | |
802 | size_type s2 = stride (c2 ()); | |
803 | result_type t = result_type (0); | |
804 | if (s1 == 1 && s2 == 1) { | |
805 | for (size_type j = 0; j < size; ++ j) | |
806 | t += data1 [j] * data2 [j]; | |
807 | } else if (s2 == 1) { | |
808 | for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1) | |
809 | t += data1 [j1] * data2 [j]; | |
810 | } else if (s1 == 1) { | |
811 | for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2) | |
812 | t += data1 [j] * data2 [j2]; | |
813 | } else { | |
814 | for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2) | |
815 | t += data1 [j1] * data2 [j2]; | |
816 | } | |
817 | return t; | |
818 | #elif defined(BOOST_UBLAS_HAVE_BINDINGS) | |
819 | return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ()); | |
820 | #else | |
821 | return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i)); | |
822 | #endif | |
823 | } | |
824 | template<class E1, class E2> | |
825 | static BOOST_UBLAS_INLINE | |
826 | result_type apply (const matrix_expression<E1> &e1, | |
827 | const vector_expression<E2> &e2, | |
828 | size_type i) { | |
829 | size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ()); | |
830 | result_type t = result_type (0); | |
831 | #ifndef BOOST_UBLAS_USE_DUFF_DEVICE | |
832 | for (size_type j = 0; j < size; ++ j) | |
833 | t += e1 () (i, j) * e2 () (j); | |
834 | #else | |
835 | size_type j (0); | |
836 | DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j)); | |
837 | #endif | |
838 | return t; | |
839 | } | |
840 | // Dense case | |
841 | template<class I1, class I2> | |
842 | static BOOST_UBLAS_INLINE | |
843 | result_type apply (difference_type size, I1 it1, I2 it2) { | |
844 | result_type t = result_type (0); | |
845 | #ifndef BOOST_UBLAS_USE_DUFF_DEVICE | |
846 | while (-- size >= 0) | |
847 | t += *it1 * *it2, ++ it1, ++ it2; | |
848 | #else | |
849 | DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); | |
850 | #endif | |
851 | return t; | |
852 | } | |
853 | // Packed case | |
854 | template<class I1, class I2> | |
855 | static BOOST_UBLAS_INLINE | |
856 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { | |
857 | result_type t = result_type (0); | |
858 | difference_type it1_size (it1_end - it1); | |
859 | difference_type it2_size (it2_end - it2); | |
860 | difference_type diff (0); | |
861 | if (it1_size > 0 && it2_size > 0) | |
862 | diff = it2.index () - it1.index2 (); | |
863 | if (diff != 0) { | |
864 | difference_type size = (std::min) (diff, it1_size); | |
865 | if (size > 0) { | |
866 | it1 += size; | |
867 | it1_size -= size; | |
868 | diff -= size; | |
869 | } | |
870 | size = (std::min) (- diff, it2_size); | |
871 | if (size > 0) { | |
872 | it2 += size; | |
873 | it2_size -= size; | |
874 | diff += size; | |
875 | } | |
876 | } | |
877 | difference_type size ((std::min) (it1_size, it2_size)); | |
878 | while (-- size >= 0) | |
879 | t += *it1 * *it2, ++ it1, ++ it2; | |
880 | return t; | |
881 | } | |
882 | // Sparse case | |
883 | template<class I1, class I2> | |
884 | static BOOST_UBLAS_INLINE | |
885 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, | |
886 | sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) { | |
887 | result_type t = result_type (0); | |
888 | if (it1 != it1_end && it2 != it2_end) { | |
889 | size_type it1_index = it1.index2 (), it2_index = it2.index (); | |
890 | while (true) { | |
891 | difference_type compare = it1_index - it2_index; | |
892 | if (compare == 0) { | |
893 | t += *it1 * *it2, ++ it1, ++ it2; | |
894 | if (it1 != it1_end && it2 != it2_end) { | |
895 | it1_index = it1.index2 (); | |
896 | it2_index = it2.index (); | |
897 | } else | |
898 | break; | |
899 | } else if (compare < 0) { | |
900 | increment (it1, it1_end, - compare); | |
901 | if (it1 != it1_end) | |
902 | it1_index = it1.index2 (); | |
903 | else | |
904 | break; | |
905 | } else if (compare > 0) { | |
906 | increment (it2, it2_end, compare); | |
907 | if (it2 != it2_end) | |
908 | it2_index = it2.index (); | |
909 | else | |
910 | break; | |
911 | } | |
912 | } | |
913 | } | |
914 | return t; | |
915 | } | |
916 | // Sparse packed case | |
917 | template<class I1, class I2> | |
918 | static BOOST_UBLAS_INLINE | |
919 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */, | |
920 | sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) { | |
921 | result_type t = result_type (0); | |
922 | while (it1 != it1_end) { | |
923 | t += *it1 * it2 () (it1.index2 ()); | |
924 | ++ it1; | |
925 | } | |
926 | return t; | |
927 | } | |
928 | // Packed sparse case | |
929 | template<class I1, class I2> | |
930 | static BOOST_UBLAS_INLINE | |
931 | result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end, | |
932 | packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) { | |
933 | result_type t = result_type (0); | |
934 | while (it2 != it2_end) { | |
935 | t += it1 () (it1.index1 (), it2.index ()) * *it2; | |
936 | ++ it2; | |
937 | } | |
938 | return t; | |
939 | } | |
940 | // Another dispatcher | |
941 | template<class I1, class I2> | |
942 | static BOOST_UBLAS_INLINE | |
943 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, | |
944 | sparse_bidirectional_iterator_tag) { | |
945 | typedef typename I1::iterator_category iterator1_category; | |
946 | typedef typename I2::iterator_category iterator2_category; | |
947 | return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ()); | |
948 | } | |
949 | }; | |
950 | ||
951 | template<class M1, class M2, class TV> | |
952 | struct matrix_vector_prod2: | |
953 | public matrix_vector_binary_functor<M1, M2, TV> { | |
954 | typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type; | |
955 | typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type; | |
956 | typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type; | |
957 | typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type; | |
958 | ||
959 | template<class C1, class C2> | |
960 | static BOOST_UBLAS_INLINE | |
961 | result_type apply (const vector_container<C1> &c1, | |
962 | const matrix_container<C2> &c2, | |
963 | size_type i) { | |
964 | #ifdef BOOST_UBLAS_USE_SIMD | |
965 | using namespace raw; | |
966 | size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ()); | |
967 | const typename M1::value_type *data1 = data_const (c1 ()); | |
968 | const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ()); | |
969 | size_type s1 = stride (c1 ()); | |
970 | size_type s2 = stride1 (c2 ()); | |
971 | result_type t = result_type (0); | |
972 | if (s1 == 1 && s2 == 1) { | |
973 | for (size_type j = 0; j < size; ++ j) | |
974 | t += data1 [j] * data2 [j]; | |
975 | } else if (s2 == 1) { | |
976 | for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1) | |
977 | t += data1 [j1] * data2 [j]; | |
978 | } else if (s1 == 1) { | |
979 | for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2) | |
980 | t += data1 [j] * data2 [j2]; | |
981 | } else { | |
982 | for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2) | |
983 | t += data1 [j1] * data2 [j2]; | |
984 | } | |
985 | return t; | |
986 | #elif defined(BOOST_UBLAS_HAVE_BINDINGS) | |
987 | return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i)); | |
988 | #else | |
989 | return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i)); | |
990 | #endif | |
991 | } | |
992 | template<class E1, class E2> | |
993 | static BOOST_UBLAS_INLINE | |
994 | result_type apply (const vector_expression<E1> &e1, | |
995 | const matrix_expression<E2> &e2, | |
996 | size_type i) { | |
997 | size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ()); | |
998 | result_type t = result_type (0); | |
999 | #ifndef BOOST_UBLAS_USE_DUFF_DEVICE | |
1000 | for (size_type j = 0; j < size; ++ j) | |
1001 | t += e1 () (j) * e2 () (j, i); | |
1002 | #else | |
1003 | size_type j (0); | |
1004 | DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j)); | |
1005 | #endif | |
1006 | return t; | |
1007 | } | |
1008 | // Dense case | |
1009 | template<class I1, class I2> | |
1010 | static BOOST_UBLAS_INLINE | |
1011 | result_type apply (difference_type size, I1 it1, I2 it2) { | |
1012 | result_type t = result_type (0); | |
1013 | #ifndef BOOST_UBLAS_USE_DUFF_DEVICE | |
1014 | while (-- size >= 0) | |
1015 | t += *it1 * *it2, ++ it1, ++ it2; | |
1016 | #else | |
1017 | DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); | |
1018 | #endif | |
1019 | return t; | |
1020 | } | |
1021 | // Packed case | |
1022 | template<class I1, class I2> | |
1023 | static BOOST_UBLAS_INLINE | |
1024 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { | |
1025 | result_type t = result_type (0); | |
1026 | difference_type it1_size (it1_end - it1); | |
1027 | difference_type it2_size (it2_end - it2); | |
1028 | difference_type diff (0); | |
1029 | if (it1_size > 0 && it2_size > 0) | |
1030 | diff = it2.index1 () - it1.index (); | |
1031 | if (diff != 0) { | |
1032 | difference_type size = (std::min) (diff, it1_size); | |
1033 | if (size > 0) { | |
1034 | it1 += size; | |
1035 | it1_size -= size; | |
1036 | diff -= size; | |
1037 | } | |
1038 | size = (std::min) (- diff, it2_size); | |
1039 | if (size > 0) { | |
1040 | it2 += size; | |
1041 | it2_size -= size; | |
1042 | diff += size; | |
1043 | } | |
1044 | } | |
1045 | difference_type size ((std::min) (it1_size, it2_size)); | |
1046 | while (-- size >= 0) | |
1047 | t += *it1 * *it2, ++ it1, ++ it2; | |
1048 | return t; | |
1049 | } | |
1050 | // Sparse case | |
1051 | template<class I1, class I2> | |
1052 | static BOOST_UBLAS_INLINE | |
1053 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, | |
1054 | sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) { | |
1055 | result_type t = result_type (0); | |
1056 | if (it1 != it1_end && it2 != it2_end) { | |
1057 | size_type it1_index = it1.index (), it2_index = it2.index1 (); | |
1058 | while (true) { | |
1059 | difference_type compare = it1_index - it2_index; | |
1060 | if (compare == 0) { | |
1061 | t += *it1 * *it2, ++ it1, ++ it2; | |
1062 | if (it1 != it1_end && it2 != it2_end) { | |
1063 | it1_index = it1.index (); | |
1064 | it2_index = it2.index1 (); | |
1065 | } else | |
1066 | break; | |
1067 | } else if (compare < 0) { | |
1068 | increment (it1, it1_end, - compare); | |
1069 | if (it1 != it1_end) | |
1070 | it1_index = it1.index (); | |
1071 | else | |
1072 | break; | |
1073 | } else if (compare > 0) { | |
1074 | increment (it2, it2_end, compare); | |
1075 | if (it2 != it2_end) | |
1076 | it2_index = it2.index1 (); | |
1077 | else | |
1078 | break; | |
1079 | } | |
1080 | } | |
1081 | } | |
1082 | return t; | |
1083 | } | |
1084 | // Packed sparse case | |
1085 | template<class I1, class I2> | |
1086 | static BOOST_UBLAS_INLINE | |
1087 | result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end, | |
1088 | packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) { | |
1089 | result_type t = result_type (0); | |
1090 | while (it2 != it2_end) { | |
1091 | t += it1 () (it2.index1 ()) * *it2; | |
1092 | ++ it2; | |
1093 | } | |
1094 | return t; | |
1095 | } | |
1096 | // Sparse packed case | |
1097 | template<class I1, class I2> | |
1098 | static BOOST_UBLAS_INLINE | |
1099 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */, | |
1100 | sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) { | |
1101 | result_type t = result_type (0); | |
1102 | while (it1 != it1_end) { | |
1103 | t += *it1 * it2 () (it1.index (), it2.index2 ()); | |
1104 | ++ it1; | |
1105 | } | |
1106 | return t; | |
1107 | } | |
1108 | // Another dispatcher | |
1109 | template<class I1, class I2> | |
1110 | static BOOST_UBLAS_INLINE | |
1111 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, | |
1112 | sparse_bidirectional_iterator_tag) { | |
1113 | typedef typename I1::iterator_category iterator1_category; | |
1114 | typedef typename I2::iterator_category iterator2_category; | |
1115 | return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ()); | |
1116 | } | |
1117 | }; | |
1118 | ||
1119 | // Binary returning matrix | |
1120 | template<class M1, class M2, class TV> | |
1121 | struct matrix_matrix_binary_functor { | |
1122 | typedef typename M1::size_type size_type; | |
1123 | typedef typename M1::difference_type difference_type; | |
1124 | typedef TV value_type; | |
1125 | typedef TV result_type; | |
1126 | }; | |
1127 | ||
1128 | template<class M1, class M2, class TV> | |
1129 | struct matrix_matrix_prod: | |
1130 | public matrix_matrix_binary_functor<M1, M2, TV> { | |
1131 | typedef typename matrix_matrix_binary_functor<M1, M2, TV>::size_type size_type; | |
1132 | typedef typename matrix_matrix_binary_functor<M1, M2, TV>::difference_type difference_type; | |
1133 | typedef typename matrix_matrix_binary_functor<M1, M2, TV>::value_type value_type; | |
1134 | typedef typename matrix_matrix_binary_functor<M1, M2, TV>::result_type result_type; | |
1135 | ||
1136 | template<class C1, class C2> | |
1137 | static BOOST_UBLAS_INLINE | |
1138 | result_type apply (const matrix_container<C1> &c1, | |
1139 | const matrix_container<C2> &c2, | |
1140 | size_type i, size_type j) { | |
1141 | #ifdef BOOST_UBLAS_USE_SIMD | |
1142 | using namespace raw; | |
1143 | size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ()); | |
1144 | const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ()); | |
1145 | const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ()); | |
1146 | size_type s1 = stride2 (c1 ()); | |
1147 | size_type s2 = stride1 (c2 ()); | |
1148 | result_type t = result_type (0); | |
1149 | if (s1 == 1 && s2 == 1) { | |
1150 | for (size_type k = 0; k < size; ++ k) | |
1151 | t += data1 [k] * data2 [k]; | |
1152 | } else if (s2 == 1) { | |
1153 | for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1) | |
1154 | t += data1 [k1] * data2 [k]; | |
1155 | } else if (s1 == 1) { | |
1156 | for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2) | |
1157 | t += data1 [k] * data2 [k2]; | |
1158 | } else { | |
1159 | for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2) | |
1160 | t += data1 [k1] * data2 [k2]; | |
1161 | } | |
1162 | return t; | |
1163 | #elif defined(BOOST_UBLAS_HAVE_BINDINGS) | |
1164 | return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j)); | |
1165 | #else | |
1166 | boost::ignore_unused(j); | |
1167 | return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i)); | |
1168 | #endif | |
1169 | } | |
1170 | template<class E1, class E2> | |
1171 | static BOOST_UBLAS_INLINE | |
1172 | result_type apply (const matrix_expression<E1> &e1, | |
1173 | const matrix_expression<E2> &e2, | |
1174 | size_type i, size_type j) { | |
1175 | size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()); | |
1176 | result_type t = result_type (0); | |
1177 | #ifndef BOOST_UBLAS_USE_DUFF_DEVICE | |
1178 | for (size_type k = 0; k < size; ++ k) | |
1179 | t += e1 () (i, k) * e2 () (k, j); | |
1180 | #else | |
1181 | size_type k (0); | |
1182 | DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k)); | |
1183 | #endif | |
1184 | return t; | |
1185 | } | |
1186 | // Dense case | |
1187 | template<class I1, class I2> | |
1188 | static BOOST_UBLAS_INLINE | |
1189 | result_type apply (difference_type size, I1 it1, I2 it2) { | |
1190 | result_type t = result_type (0); | |
1191 | #ifndef BOOST_UBLAS_USE_DUFF_DEVICE | |
1192 | while (-- size >= 0) | |
1193 | t += *it1 * *it2, ++ it1, ++ it2; | |
1194 | #else | |
1195 | DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); | |
1196 | #endif | |
1197 | return t; | |
1198 | } | |
1199 | // Packed case | |
1200 | template<class I1, class I2> | |
1201 | static BOOST_UBLAS_INLINE | |
1202 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) { | |
1203 | result_type t = result_type (0); | |
1204 | difference_type it1_size (it1_end - it1); | |
1205 | difference_type it2_size (it2_end - it2); | |
1206 | difference_type diff (0); | |
1207 | if (it1_size > 0 && it2_size > 0) | |
1208 | diff = it2.index1 () - it1.index2 (); | |
1209 | if (diff != 0) { | |
1210 | difference_type size = (std::min) (diff, it1_size); | |
1211 | if (size > 0) { | |
1212 | it1 += size; | |
1213 | it1_size -= size; | |
1214 | diff -= size; | |
1215 | } | |
1216 | size = (std::min) (- diff, it2_size); | |
1217 | if (size > 0) { | |
1218 | it2 += size; | |
1219 | it2_size -= size; | |
1220 | diff += size; | |
1221 | } | |
1222 | } | |
1223 | difference_type size ((std::min) (it1_size, it2_size)); | |
1224 | while (-- size >= 0) | |
1225 | t += *it1 * *it2, ++ it1, ++ it2; | |
1226 | return t; | |
1227 | } | |
1228 | // Sparse case | |
1229 | template<class I1, class I2> | |
1230 | static BOOST_UBLAS_INLINE | |
1231 | result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { | |
1232 | result_type t = result_type (0); | |
1233 | if (it1 != it1_end && it2 != it2_end) { | |
1234 | size_type it1_index = it1.index2 (), it2_index = it2.index1 (); | |
1235 | while (true) { | |
1236 | difference_type compare = difference_type (it1_index - it2_index); | |
1237 | if (compare == 0) { | |
1238 | t += *it1 * *it2, ++ it1, ++ it2; | |
1239 | if (it1 != it1_end && it2 != it2_end) { | |
1240 | it1_index = it1.index2 (); | |
1241 | it2_index = it2.index1 (); | |
1242 | } else | |
1243 | break; | |
1244 | } else if (compare < 0) { | |
1245 | increment (it1, it1_end, - compare); | |
1246 | if (it1 != it1_end) | |
1247 | it1_index = it1.index2 (); | |
1248 | else | |
1249 | break; | |
1250 | } else if (compare > 0) { | |
1251 | increment (it2, it2_end, compare); | |
1252 | if (it2 != it2_end) | |
1253 | it2_index = it2.index1 (); | |
1254 | else | |
1255 | break; | |
1256 | } | |
1257 | } | |
1258 | } | |
1259 | return t; | |
1260 | } | |
1261 | }; | |
1262 | ||
1263 | // Unary returning scalar norm | |
1264 | template<class M> | |
1265 | struct matrix_scalar_real_unary_functor { | |
1266 | typedef typename M::value_type value_type; | |
1267 | typedef typename type_traits<value_type>::real_type real_type; | |
1268 | typedef real_type result_type; | |
1269 | }; | |
1270 | ||
1271 | template<class M> | |
1272 | struct matrix_norm_1: | |
1273 | public matrix_scalar_real_unary_functor<M> { | |
1274 | typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; | |
1275 | typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; | |
1276 | typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; | |
1277 | ||
1278 | template<class E> | |
1279 | static BOOST_UBLAS_INLINE | |
1280 | result_type apply (const matrix_expression<E> &e) { | |
1281 | real_type t = real_type (); | |
1282 | typedef typename E::size_type matrix_size_type; | |
1283 | matrix_size_type size2 (e ().size2 ()); | |
1284 | for (matrix_size_type j = 0; j < size2; ++ j) { | |
1285 | real_type u = real_type (); | |
1286 | matrix_size_type size1 (e ().size1 ()); | |
1287 | for (matrix_size_type i = 0; i < size1; ++ i) { | |
1288 | real_type v (type_traits<value_type>::norm_1 (e () (i, j))); | |
1289 | u += v; | |
1290 | } | |
1291 | if (u > t) | |
1292 | t = u; | |
1293 | } | |
1294 | return t; | |
1295 | } | |
1296 | }; | |
1297 | ||
1298 | template<class M> | |
1299 | struct matrix_norm_frobenius: | |
1300 | public matrix_scalar_real_unary_functor<M> { | |
1301 | typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; | |
1302 | typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; | |
1303 | typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; | |
1304 | ||
1305 | template<class E> | |
1306 | static BOOST_UBLAS_INLINE | |
1307 | result_type apply (const matrix_expression<E> &e) { | |
1308 | real_type t = real_type (); | |
1309 | typedef typename E::size_type matrix_size_type; | |
1310 | matrix_size_type size1 (e ().size1 ()); | |
1311 | for (matrix_size_type i = 0; i < size1; ++ i) { | |
1312 | matrix_size_type size2 (e ().size2 ()); | |
1313 | for (matrix_size_type j = 0; j < size2; ++ j) { | |
1314 | real_type u (type_traits<value_type>::norm_2 (e () (i, j))); | |
1315 | t += u * u; | |
1316 | } | |
1317 | } | |
1318 | return type_traits<real_type>::type_sqrt (t); | |
1319 | } | |
1320 | }; | |
1321 | ||
1322 | template<class M> | |
1323 | struct matrix_norm_inf: | |
1324 | public matrix_scalar_real_unary_functor<M> { | |
1325 | typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; | |
1326 | typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; | |
1327 | typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; | |
1328 | ||
1329 | template<class E> | |
1330 | static BOOST_UBLAS_INLINE | |
1331 | result_type apply (const matrix_expression<E> &e) { | |
1332 | real_type t = real_type (); | |
1333 | typedef typename E::size_type matrix_size_type; | |
1334 | matrix_size_type size1 (e ().size1 ()); | |
1335 | for (matrix_size_type i = 0; i < size1; ++ i) { | |
1336 | real_type u = real_type (); | |
1337 | matrix_size_type size2 (e ().size2 ()); | |
1338 | for (matrix_size_type j = 0; j < size2; ++ j) { | |
1339 | real_type v (type_traits<value_type>::norm_inf (e () (i, j))); | |
1340 | u += v; | |
1341 | } | |
1342 | if (u > t) | |
1343 | t = u; | |
1344 | } | |
1345 | return t; | |
1346 | } | |
1347 | }; | |
1348 | ||
1349 | // forward declaration | |
1350 | template <class Z, class D> struct basic_column_major; | |
1351 | ||
1352 | // This functor defines storage layout and it's properties | |
1353 | // matrix (i,j) -> storage [i * size_i + j] | |
1354 | template <class Z, class D> | |
1355 | struct basic_row_major { | |
1356 | typedef Z size_type; | |
1357 | typedef D difference_type; | |
1358 | typedef row_major_tag orientation_category; | |
1359 | typedef basic_column_major<Z,D> transposed_layout; | |
1360 | ||
1361 | static | |
1362 | BOOST_UBLAS_INLINE | |
1363 | size_type storage_size (size_type size_i, size_type size_j) { | |
1364 | // Guard against size_type overflow | |
1365 | BOOST_UBLAS_CHECK (size_j == 0 || size_i <= (std::numeric_limits<size_type>::max) () / size_j, bad_size ()); | |
1366 | return size_i * size_j; | |
1367 | } | |
1368 | ||
1369 | // Indexing conversion to storage element | |
1370 | static | |
1371 | BOOST_UBLAS_INLINE | |
1372 | size_type element (size_type i, size_type size_i, size_type j, size_type size_j) { | |
1373 | BOOST_UBLAS_CHECK (i < size_i, bad_index ()); | |
1374 | BOOST_UBLAS_CHECK (j < size_j, bad_index ()); | |
1375 | detail::ignore_unused_variable_warning(size_i); | |
1376 | // Guard against size_type overflow | |
1377 | BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ()); | |
1378 | return i * size_j + j; | |
1379 | } | |
1380 | static | |
1381 | BOOST_UBLAS_INLINE | |
1382 | size_type address (size_type i, size_type size_i, size_type j, size_type size_j) { | |
1383 | BOOST_UBLAS_CHECK (i <= size_i, bad_index ()); | |
1384 | BOOST_UBLAS_CHECK (j <= size_j, bad_index ()); | |
1385 | // Guard against size_type overflow - address may be size_j past end of storage | |
1386 | BOOST_UBLAS_CHECK (size_j == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ()); | |
1387 | detail::ignore_unused_variable_warning(size_i); | |
1388 | return i * size_j + j; | |
1389 | } | |
1390 | ||
1391 | // Storage element to index conversion | |
1392 | static | |
1393 | BOOST_UBLAS_INLINE | |
1394 | difference_type distance_i (difference_type k, size_type /* size_i */, size_type size_j) { | |
1395 | return size_j != 0 ? k / size_j : 0; | |
1396 | } | |
1397 | static | |
1398 | BOOST_UBLAS_INLINE | |
1399 | difference_type distance_j (difference_type k, size_type /* size_i */, size_type /* size_j */) { | |
1400 | return k; | |
1401 | } | |
1402 | static | |
1403 | BOOST_UBLAS_INLINE | |
1404 | size_type index_i (difference_type k, size_type /* size_i */, size_type size_j) { | |
1405 | return size_j != 0 ? k / size_j : 0; | |
1406 | } | |
1407 | static | |
1408 | BOOST_UBLAS_INLINE | |
1409 | size_type index_j (difference_type k, size_type /* size_i */, size_type size_j) { | |
1410 | return size_j != 0 ? k % size_j : 0; | |
1411 | } | |
1412 | static | |
1413 | BOOST_UBLAS_INLINE | |
1414 | bool fast_i () { | |
1415 | return false; | |
1416 | } | |
1417 | static | |
1418 | BOOST_UBLAS_INLINE | |
1419 | bool fast_j () { | |
1420 | return true; | |
1421 | } | |
1422 | ||
1423 | // Iterating storage elements | |
1424 | template<class I> | |
1425 | static | |
1426 | BOOST_UBLAS_INLINE | |
1427 | void increment_i (I &it, size_type /* size_i */, size_type size_j) { | |
1428 | it += size_j; | |
1429 | } | |
1430 | template<class I> | |
1431 | static | |
1432 | BOOST_UBLAS_INLINE | |
1433 | void increment_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) { | |
1434 | it += n * size_j; | |
1435 | } | |
1436 | template<class I> | |
1437 | static | |
1438 | BOOST_UBLAS_INLINE | |
1439 | void decrement_i (I &it, size_type /* size_i */, size_type size_j) { | |
1440 | it -= size_j; | |
1441 | } | |
1442 | template<class I> | |
1443 | static | |
1444 | BOOST_UBLAS_INLINE | |
1445 | void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) { | |
1446 | it -= n * size_j; | |
1447 | } | |
1448 | template<class I> | |
1449 | static | |
1450 | BOOST_UBLAS_INLINE | |
1451 | void increment_j (I &it, size_type /* size_i */, size_type /* size_j */) { | |
1452 | ++ it; | |
1453 | } | |
1454 | template<class I> | |
1455 | static | |
1456 | BOOST_UBLAS_INLINE | |
1457 | void increment_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { | |
1458 | it += n; | |
1459 | } | |
1460 | template<class I> | |
1461 | static | |
1462 | BOOST_UBLAS_INLINE | |
1463 | void decrement_j (I &it, size_type /* size_i */, size_type /* size_j */) { | |
1464 | -- it; | |
1465 | } | |
1466 | template<class I> | |
1467 | static | |
1468 | BOOST_UBLAS_INLINE | |
1469 | void decrement_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { | |
1470 | it -= n; | |
1471 | } | |
1472 | ||
1473 | // Triangular access | |
1474 | static | |
1475 | BOOST_UBLAS_INLINE | |
1476 | size_type triangular_size (size_type size_i, size_type size_j) { | |
1477 | size_type size = (std::max) (size_i, size_j); | |
1478 | // Guard against size_type overflow - simplified | |
1479 | BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ()); | |
1480 | return ((size + 1) * size) / 2; | |
1481 | } | |
1482 | static | |
1483 | BOOST_UBLAS_INLINE | |
1484 | size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) { | |
1485 | BOOST_UBLAS_CHECK (i < size_i, bad_index ()); | |
1486 | BOOST_UBLAS_CHECK (j < size_j, bad_index ()); | |
1487 | BOOST_UBLAS_CHECK (i >= j, bad_index ()); | |
1488 | detail::ignore_unused_variable_warning(size_i); | |
1489 | detail::ignore_unused_variable_warning(size_j); | |
1490 | // FIXME size_type overflow | |
1491 | // sigma_i (i + 1) = (i + 1) * i / 2 | |
1492 | // i = 0 1 2 3, sigma = 0 1 3 6 | |
1493 | return ((i + 1) * i) / 2 + j; | |
1494 | } | |
1495 | static | |
1496 | BOOST_UBLAS_INLINE | |
1497 | size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) { | |
1498 | BOOST_UBLAS_CHECK (i < size_i, bad_index ()); | |
1499 | BOOST_UBLAS_CHECK (j < size_j, bad_index ()); | |
1500 | BOOST_UBLAS_CHECK (i <= j, bad_index ()); | |
1501 | // FIXME size_type overflow | |
1502 | // sigma_i (size - i) = size * i - i * (i - 1) / 2 | |
1503 | // i = 0 1 2 3, sigma = 0 4 7 9 | |
1504 | return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i; | |
1505 | } | |
1506 | ||
1507 | // Major and minor indices | |
1508 | static | |
1509 | BOOST_UBLAS_INLINE | |
1510 | size_type index_M (size_type index1, size_type /* index2 */) { | |
1511 | return index1; | |
1512 | } | |
1513 | static | |
1514 | BOOST_UBLAS_INLINE | |
1515 | size_type index_m (size_type /* index1 */, size_type index2) { | |
1516 | return index2; | |
1517 | } | |
1518 | static | |
1519 | BOOST_UBLAS_INLINE | |
1520 | size_type size_M (size_type size_i, size_type /* size_j */) { | |
1521 | return size_i; | |
1522 | } | |
1523 | static | |
1524 | BOOST_UBLAS_INLINE | |
1525 | size_type size_m (size_type /* size_i */, size_type size_j) { | |
1526 | return size_j; | |
1527 | } | |
1528 | }; | |
1529 | ||
1530 | // This functor defines storage layout and it's properties | |
1531 | // matrix (i,j) -> storage [i + j * size_i] | |
1532 | template <class Z, class D> | |
1533 | struct basic_column_major { | |
1534 | typedef Z size_type; | |
1535 | typedef D difference_type; | |
1536 | typedef column_major_tag orientation_category; | |
1537 | typedef basic_row_major<Z,D> transposed_layout; | |
1538 | ||
1539 | static | |
1540 | BOOST_UBLAS_INLINE | |
1541 | size_type storage_size (size_type size_i, size_type size_j) { | |
1542 | // Guard against size_type overflow | |
1543 | BOOST_UBLAS_CHECK (size_i == 0 || size_j <= (std::numeric_limits<size_type>::max) () / size_i, bad_size ()); | |
1544 | return size_i * size_j; | |
1545 | } | |
1546 | ||
1547 | // Indexing conversion to storage element | |
1548 | static | |
1549 | BOOST_UBLAS_INLINE | |
1550 | size_type element (size_type i, size_type size_i, size_type j, size_type size_j) { | |
1551 | BOOST_UBLAS_CHECK (i < size_i, bad_index ()); | |
1552 | BOOST_UBLAS_CHECK (j < size_j, bad_index ()); | |
1553 | detail::ignore_unused_variable_warning(size_j); | |
1554 | // Guard against size_type overflow | |
1555 | BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ()); | |
1556 | return i + j * size_i; | |
1557 | } | |
1558 | static | |
1559 | BOOST_UBLAS_INLINE | |
1560 | size_type address (size_type i, size_type size_i, size_type j, size_type size_j) { | |
1561 | BOOST_UBLAS_CHECK (i <= size_i, bad_index ()); | |
1562 | BOOST_UBLAS_CHECK (j <= size_j, bad_index ()); | |
1563 | detail::ignore_unused_variable_warning(size_j); | |
1564 | // Guard against size_type overflow - address may be size_i past end of storage | |
1565 | BOOST_UBLAS_CHECK (size_i == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ()); | |
1566 | return i + j * size_i; | |
1567 | } | |
1568 | ||
1569 | // Storage element to index conversion | |
1570 | static | |
1571 | BOOST_UBLAS_INLINE | |
1572 | difference_type distance_i (difference_type k, size_type /* size_i */, size_type /* size_j */) { | |
1573 | return k; | |
1574 | } | |
1575 | static | |
1576 | BOOST_UBLAS_INLINE | |
1577 | difference_type distance_j (difference_type k, size_type size_i, size_type /* size_j */) { | |
1578 | return size_i != 0 ? k / size_i : 0; | |
1579 | } | |
1580 | static | |
1581 | BOOST_UBLAS_INLINE | |
1582 | size_type index_i (difference_type k, size_type size_i, size_type /* size_j */) { | |
1583 | return size_i != 0 ? k % size_i : 0; | |
1584 | } | |
1585 | static | |
1586 | BOOST_UBLAS_INLINE | |
1587 | size_type index_j (difference_type k, size_type size_i, size_type /* size_j */) { | |
1588 | return size_i != 0 ? k / size_i : 0; | |
1589 | } | |
1590 | static | |
1591 | BOOST_UBLAS_INLINE | |
1592 | bool fast_i () { | |
1593 | return true; | |
1594 | } | |
1595 | static | |
1596 | BOOST_UBLAS_INLINE | |
1597 | bool fast_j () { | |
1598 | return false; | |
1599 | } | |
1600 | ||
1601 | // Iterating | |
1602 | template<class I> | |
1603 | static | |
1604 | BOOST_UBLAS_INLINE | |
1605 | void increment_i (I &it, size_type /* size_i */, size_type /* size_j */) { | |
1606 | ++ it; | |
1607 | } | |
1608 | template<class I> | |
1609 | static | |
1610 | BOOST_UBLAS_INLINE | |
1611 | void increment_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { | |
1612 | it += n; | |
1613 | } | |
1614 | template<class I> | |
1615 | static | |
1616 | BOOST_UBLAS_INLINE | |
1617 | void decrement_i (I &it, size_type /* size_i */, size_type /* size_j */) { | |
1618 | -- it; | |
1619 | } | |
1620 | template<class I> | |
1621 | static | |
1622 | BOOST_UBLAS_INLINE | |
1623 | void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { | |
1624 | it -= n; | |
1625 | } | |
1626 | template<class I> | |
1627 | static | |
1628 | BOOST_UBLAS_INLINE | |
1629 | void increment_j (I &it, size_type size_i, size_type /* size_j */) { | |
1630 | it += size_i; | |
1631 | } | |
1632 | template<class I> | |
1633 | static | |
1634 | BOOST_UBLAS_INLINE | |
1635 | void increment_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) { | |
1636 | it += n * size_i; | |
1637 | } | |
1638 | template<class I> | |
1639 | static | |
1640 | BOOST_UBLAS_INLINE | |
1641 | void decrement_j (I &it, size_type size_i, size_type /* size_j */) { | |
1642 | it -= size_i; | |
1643 | } | |
1644 | template<class I> | |
1645 | static | |
1646 | BOOST_UBLAS_INLINE | |
1647 | void decrement_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) { | |
1648 | it -= n* size_i; | |
1649 | } | |
1650 | ||
1651 | // Triangular access | |
1652 | static | |
1653 | BOOST_UBLAS_INLINE | |
1654 | size_type triangular_size (size_type size_i, size_type size_j) { | |
1655 | size_type size = (std::max) (size_i, size_j); | |
1656 | // Guard against size_type overflow - simplified | |
1657 | BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ()); | |
1658 | return ((size + 1) * size) / 2; | |
1659 | } | |
1660 | static | |
1661 | BOOST_UBLAS_INLINE | |
1662 | size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) { | |
1663 | BOOST_UBLAS_CHECK (i < size_i, bad_index ()); | |
1664 | BOOST_UBLAS_CHECK (j < size_j, bad_index ()); | |
1665 | BOOST_UBLAS_CHECK (i >= j, bad_index ()); | |
1666 | // FIXME size_type overflow | |
1667 | // sigma_j (size - j) = size * j - j * (j - 1) / 2 | |
1668 | // j = 0 1 2 3, sigma = 0 4 7 9 | |
1669 | return i - j + (j * (2 * (std::max) (size_i, size_j) - j + 1)) / 2; | |
1670 | } | |
1671 | static | |
1672 | BOOST_UBLAS_INLINE | |
1673 | size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) { | |
1674 | BOOST_UBLAS_CHECK (i < size_i, bad_index ()); | |
1675 | BOOST_UBLAS_CHECK (j < size_j, bad_index ()); | |
1676 | BOOST_UBLAS_CHECK (i <= j, bad_index ()); | |
1677 | // FIXME size_type overflow | |
1678 | // sigma_j (j + 1) = (j + 1) * j / 2 | |
1679 | // j = 0 1 2 3, sigma = 0 1 3 6 | |
1680 | return i + ((j + 1) * j) / 2; | |
1681 | } | |
1682 | ||
1683 | // Major and minor indices | |
1684 | static | |
1685 | BOOST_UBLAS_INLINE | |
1686 | size_type index_M (size_type /* index1 */, size_type index2) { | |
1687 | return index2; | |
1688 | } | |
1689 | static | |
1690 | BOOST_UBLAS_INLINE | |
1691 | size_type index_m (size_type index1, size_type /* index2 */) { | |
1692 | return index1; | |
1693 | } | |
1694 | static | |
1695 | BOOST_UBLAS_INLINE | |
1696 | size_type size_M (size_type /* size_i */, size_type size_j) { | |
1697 | return size_j; | |
1698 | } | |
1699 | static | |
1700 | BOOST_UBLAS_INLINE | |
1701 | size_type size_m (size_type size_i, size_type /* size_j */) { | |
1702 | return size_i; | |
1703 | } | |
1704 | }; | |
1705 | ||
1706 | ||
1707 | template <class Z> | |
1708 | struct basic_full { | |
1709 | typedef Z size_type; | |
1710 | ||
1711 | template<class L> | |
1712 | static | |
1713 | BOOST_UBLAS_INLINE | |
1714 | size_type packed_size (L, size_type size_i, size_type size_j) { | |
1715 | return L::storage_size (size_i, size_j); | |
1716 | } | |
1717 | ||
1718 | static | |
1719 | BOOST_UBLAS_INLINE | |
1720 | bool zero (size_type /* i */, size_type /* j */) { | |
1721 | return false; | |
1722 | } | |
1723 | static | |
1724 | BOOST_UBLAS_INLINE | |
1725 | bool one (size_type /* i */, size_type /* j */) { | |
1726 | return false; | |
1727 | } | |
1728 | static | |
1729 | BOOST_UBLAS_INLINE | |
1730 | bool other (size_type /* i */, size_type /* j */) { | |
1731 | return true; | |
1732 | } | |
1733 | // FIXME: this should not be used at all | |
1734 | static | |
1735 | BOOST_UBLAS_INLINE | |
1736 | size_type restrict1 (size_type i, size_type /* j */) { | |
1737 | return i; | |
1738 | } | |
1739 | static | |
1740 | BOOST_UBLAS_INLINE | |
1741 | size_type restrict2 (size_type /* i */, size_type j) { | |
1742 | return j; | |
1743 | } | |
1744 | static | |
1745 | BOOST_UBLAS_INLINE | |
1746 | size_type mutable_restrict1 (size_type i, size_type /* j */) { | |
1747 | return i; | |
1748 | } | |
1749 | static | |
1750 | BOOST_UBLAS_INLINE | |
1751 | size_type mutable_restrict2 (size_type /* i */, size_type j) { | |
1752 | return j; | |
1753 | } | |
1754 | }; | |
1755 | ||
1756 | namespace detail { | |
1757 | template < class L > | |
1758 | struct transposed_structure { | |
1759 | typedef typename L::size_type size_type; | |
1760 | ||
1761 | template<class LAYOUT> | |
1762 | static | |
1763 | BOOST_UBLAS_INLINE | |
1764 | size_type packed_size (LAYOUT l, size_type size_i, size_type size_j) { | |
1765 | return L::packed_size(l, size_j, size_i); | |
1766 | } | |
1767 | ||
1768 | static | |
1769 | BOOST_UBLAS_INLINE | |
1770 | bool zero (size_type i, size_type j) { | |
1771 | return L::zero(j, i); | |
1772 | } | |
1773 | static | |
1774 | BOOST_UBLAS_INLINE | |
1775 | bool one (size_type i, size_type j) { | |
1776 | return L::one(j, i); | |
1777 | } | |
1778 | static | |
1779 | BOOST_UBLAS_INLINE | |
1780 | bool other (size_type i, size_type j) { | |
1781 | return L::other(j, i); | |
1782 | } | |
1783 | template<class LAYOUT> | |
1784 | static | |
1785 | BOOST_UBLAS_INLINE | |
1786 | size_type element (LAYOUT /* l */, size_type i, size_type size_i, size_type j, size_type size_j) { | |
1787 | return L::element(typename LAYOUT::transposed_layout(), j, size_j, i, size_i); | |
1788 | } | |
1789 | ||
1790 | static | |
1791 | BOOST_UBLAS_INLINE | |
1792 | size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { | |
1793 | return L::restrict2(j, i, size2, size1); | |
1794 | } | |
1795 | static | |
1796 | BOOST_UBLAS_INLINE | |
1797 | size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) { | |
1798 | return L::restrict1(j, i, size2, size1); | |
1799 | } | |
1800 | static | |
1801 | BOOST_UBLAS_INLINE | |
1802 | size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type size2) { | |
1803 | return L::mutable_restrict2(j, i, size2, size1); | |
1804 | } | |
1805 | static | |
1806 | BOOST_UBLAS_INLINE | |
1807 | size_type mutable_restrict2 (size_type i, size_type j, size_type size1, size_type size2) { | |
1808 | return L::mutable_restrict1(j, i, size2, size1); | |
1809 | } | |
1810 | ||
1811 | static | |
1812 | BOOST_UBLAS_INLINE | |
1813 | size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { | |
1814 | return L::global_restrict2(index2, size2, index1, size1); | |
1815 | } | |
1816 | static | |
1817 | BOOST_UBLAS_INLINE | |
1818 | size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { | |
1819 | return L::global_restrict1(index2, size2, index1, size1); | |
1820 | } | |
1821 | static | |
1822 | BOOST_UBLAS_INLINE | |
1823 | size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { | |
1824 | return L::global_mutable_restrict2(index2, size2, index1, size1); | |
1825 | } | |
1826 | static | |
1827 | BOOST_UBLAS_INLINE | |
1828 | size_type global_mutable_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { | |
1829 | return L::global_mutable_restrict1(index2, size2, index1, size1); | |
1830 | } | |
1831 | }; | |
1832 | } | |
1833 | ||
1834 | template <class Z> | |
1835 | struct basic_lower { | |
1836 | typedef Z size_type; | |
1837 | typedef lower_tag triangular_type; | |
1838 | ||
1839 | template<class L> | |
1840 | static | |
1841 | BOOST_UBLAS_INLINE | |
1842 | size_type packed_size (L, size_type size_i, size_type size_j) { | |
1843 | return L::triangular_size (size_i, size_j); | |
1844 | } | |
1845 | ||
1846 | static | |
1847 | BOOST_UBLAS_INLINE | |
1848 | bool zero (size_type i, size_type j) { | |
1849 | return j > i; | |
1850 | } | |
1851 | static | |
1852 | BOOST_UBLAS_INLINE | |
1853 | bool one (size_type /* i */, size_type /* j */) { | |
1854 | return false; | |
1855 | } | |
1856 | static | |
1857 | BOOST_UBLAS_INLINE | |
1858 | bool other (size_type i, size_type j) { | |
1859 | return j <= i; | |
1860 | } | |
1861 | template<class L> | |
1862 | static | |
1863 | BOOST_UBLAS_INLINE | |
1864 | size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { | |
1865 | return L::lower_element (i, size_i, j, size_j); | |
1866 | } | |
1867 | ||
1868 | // return nearest valid index in column j | |
1869 | static | |
1870 | BOOST_UBLAS_INLINE | |
1871 | size_type restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { | |
1872 | return (std::max)(j, (std::min) (size1, i)); | |
1873 | } | |
1874 | // return nearest valid index in row i | |
1875 | static | |
1876 | BOOST_UBLAS_INLINE | |
1877 | size_type restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { | |
1878 | return (std::max)(size_type(0), (std::min) (i+1, j)); | |
1879 | } | |
1880 | // return nearest valid mutable index in column j | |
1881 | static | |
1882 | BOOST_UBLAS_INLINE | |
1883 | size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { | |
1884 | return (std::max)(j, (std::min) (size1, i)); | |
1885 | } | |
1886 | // return nearest valid mutable index in row i | |
1887 | static | |
1888 | BOOST_UBLAS_INLINE | |
1889 | size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { | |
1890 | return (std::max)(size_type(0), (std::min) (i+1, j)); | |
1891 | } | |
1892 | ||
1893 | // return an index between the first and (1+last) filled row | |
1894 | static | |
1895 | BOOST_UBLAS_INLINE | |
1896 | size_type global_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { | |
1897 | return (std::max)(size_type(0), (std::min)(size1, index1) ); | |
1898 | } | |
1899 | // return an index between the first and (1+last) filled column | |
1900 | static | |
1901 | BOOST_UBLAS_INLINE | |
1902 | size_type global_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { | |
1903 | return (std::max)(size_type(0), (std::min)(size2, index2) ); | |
1904 | } | |
1905 | ||
1906 | // return an index between the first and (1+last) filled mutable row | |
1907 | static | |
1908 | BOOST_UBLAS_INLINE | |
1909 | size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { | |
1910 | return (std::max)(size_type(0), (std::min)(size1, index1) ); | |
1911 | } | |
1912 | // return an index between the first and (1+last) filled mutable column | |
1913 | static | |
1914 | BOOST_UBLAS_INLINE | |
1915 | size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { | |
1916 | return (std::max)(size_type(0), (std::min)(size2, index2) ); | |
1917 | } | |
1918 | }; | |
1919 | ||
1920 | // the first row only contains a single 1. Thus it is not stored. | |
1921 | template <class Z> | |
1922 | struct basic_unit_lower : public basic_lower<Z> { | |
1923 | typedef Z size_type; | |
1924 | typedef unit_lower_tag triangular_type; | |
1925 | ||
1926 | template<class L> | |
1927 | static | |
1928 | BOOST_UBLAS_INLINE | |
1929 | size_type packed_size (L, size_type size_i, size_type size_j) { | |
1930 | // Zero size strict triangles are bad at this point | |
1931 | BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ()); | |
1932 | return L::triangular_size (size_i - 1, size_j - 1); | |
1933 | } | |
1934 | ||
1935 | static | |
1936 | BOOST_UBLAS_INLINE | |
1937 | bool one (size_type i, size_type j) { | |
1938 | return j == i; | |
1939 | } | |
1940 | static | |
1941 | BOOST_UBLAS_INLINE | |
1942 | bool other (size_type i, size_type j) { | |
1943 | return j < i; | |
1944 | } | |
1945 | template<class L> | |
1946 | static | |
1947 | BOOST_UBLAS_INLINE | |
1948 | size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { | |
1949 | // Zero size strict triangles are bad at this point | |
1950 | BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ()); | |
1951 | return L::lower_element (i-1, size_i - 1, j, size_j - 1); | |
1952 | } | |
1953 | ||
1954 | static | |
1955 | BOOST_UBLAS_INLINE | |
1956 | size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { | |
1957 | return (std::max)(j+1, (std::min) (size1, i)); | |
1958 | } | |
1959 | static | |
1960 | BOOST_UBLAS_INLINE | |
1961 | size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { | |
1962 | return (std::max)(size_type(0), (std::min) (i, j)); | |
1963 | } | |
1964 | ||
1965 | // return an index between the first and (1+last) filled mutable row | |
1966 | static | |
1967 | BOOST_UBLAS_INLINE | |
1968 | size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { | |
1969 | return (std::max)(size_type(1), (std::min)(size1, index1) ); | |
1970 | } | |
1971 | // return an index between the first and (1+last) filled mutable column | |
1972 | static | |
1973 | BOOST_UBLAS_INLINE | |
1974 | size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { | |
1975 | BOOST_UBLAS_CHECK( size2 >= 1 , external_logic() ); | |
1976 | return (std::max)(size_type(0), (std::min)(size2-1, index2) ); | |
1977 | } | |
1978 | }; | |
1979 | ||
1980 | // the first row only contains no element. Thus it is not stored. | |
1981 | template <class Z> | |
1982 | struct basic_strict_lower : public basic_unit_lower<Z> { | |
1983 | typedef Z size_type; | |
1984 | typedef strict_lower_tag triangular_type; | |
1985 | ||
1986 | template<class L> | |
1987 | static | |
1988 | BOOST_UBLAS_INLINE | |
1989 | size_type packed_size (L, size_type size_i, size_type size_j) { | |
1990 | // Zero size strict triangles are bad at this point | |
1991 | BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ()); | |
1992 | return L::triangular_size (size_i - 1, size_j - 1); | |
1993 | } | |
1994 | ||
1995 | static | |
1996 | BOOST_UBLAS_INLINE | |
1997 | bool zero (size_type i, size_type j) { | |
1998 | return j >= i; | |
1999 | } | |
2000 | static | |
2001 | BOOST_UBLAS_INLINE | |
2002 | bool one (size_type /*i*/, size_type /*j*/) { | |
2003 | return false; | |
2004 | } | |
2005 | static | |
2006 | BOOST_UBLAS_INLINE | |
2007 | bool other (size_type i, size_type j) { | |
2008 | return j < i; | |
2009 | } | |
2010 | template<class L> | |
2011 | static | |
2012 | BOOST_UBLAS_INLINE | |
2013 | size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { | |
2014 | // Zero size strict triangles are bad at this point | |
2015 | BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ()); | |
2016 | return L::lower_element (i-1, size_i - 1, j, size_j - 1); | |
2017 | } | |
2018 | ||
2019 | static | |
2020 | BOOST_UBLAS_INLINE | |
2021 | size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { | |
2022 | return basic_unit_lower<Z>::mutable_restrict1(i, j, size1, size2); | |
2023 | } | |
2024 | static | |
2025 | BOOST_UBLAS_INLINE | |
2026 | size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) { | |
2027 | return basic_unit_lower<Z>::mutable_restrict2(i, j, size1, size2); | |
2028 | } | |
2029 | ||
2030 | // return an index between the first and (1+last) filled row | |
2031 | static | |
2032 | BOOST_UBLAS_INLINE | |
2033 | size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { | |
2034 | return basic_unit_lower<Z>::global_mutable_restrict1(index1, size1, index2, size2); | |
2035 | } | |
2036 | // return an index between the first and (1+last) filled column | |
2037 | static | |
2038 | BOOST_UBLAS_INLINE | |
2039 | size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { | |
2040 | return basic_unit_lower<Z>::global_mutable_restrict2(index1, size1, index2, size2); | |
2041 | } | |
2042 | }; | |
2043 | ||
2044 | ||
2045 | template <class Z> | |
2046 | struct basic_upper : public detail::transposed_structure<basic_lower<Z> > | |
2047 | { | |
2048 | typedef upper_tag triangular_type; | |
2049 | }; | |
2050 | ||
2051 | template <class Z> | |
2052 | struct basic_unit_upper : public detail::transposed_structure<basic_unit_lower<Z> > | |
2053 | { | |
2054 | typedef unit_upper_tag triangular_type; | |
2055 | }; | |
2056 | ||
2057 | template <class Z> | |
2058 | struct basic_strict_upper : public detail::transposed_structure<basic_strict_lower<Z> > | |
2059 | { | |
2060 | typedef strict_upper_tag triangular_type; | |
2061 | }; | |
2062 | ||
2063 | ||
2064 | }}} | |
2065 | ||
2066 | #endif |