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
27namespace 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
41namespace 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 typedef typename E::size_type vector_size_type;
437 vector_size_type size (e ().size ());
438#ifndef BOOST_UBLAS_SCALED_NORM
439 real_type t = real_type ();
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 static_cast<result_type>(type_traits<real_type>::type_sqrt (t));
445#else
446 real_type scale = real_type ();
447 real_type sum_squares (1);
448 for (vector_size_type i = 0; i < size; ++ i) {
449 real_type u (type_traits<value_type>::norm_2 (e () (i)));
450 if ( real_type () /* zero */ == u ) continue;
451 if (scale < u) {
452 real_type v (scale / u);
453 sum_squares = sum_squares * v * v + real_type (1);
454 scale = u;
455 } else {
456 real_type v (u / scale);
457 sum_squares += v * v;
458 }
459 }
460 return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares));
461#endif
462 }
463 // Dense case
464 template<class D, class I>
465 static BOOST_UBLAS_INLINE
466 result_type apply (D size, I it) {
467#ifndef BOOST_UBLAS_SCALED_NORM
468 real_type t = real_type ();
469 while (-- size >= 0) {
470 real_type u (type_traits<value_type>::norm_2 (*it));
471 t += u * u;
472 ++ it;
473 }
474 return static_cast<result_type>(type_traits<real_type>::type_sqrt (t));
475#else
476 real_type scale = real_type ();
477 real_type sum_squares (1);
478 while (-- size >= 0) {
479 real_type u (type_traits<value_type>::norm_2 (*it));
480 if (scale < u) {
481 real_type v (scale / u);
482 sum_squares = sum_squares * v * v + real_type (1);
483 scale = u;
484 } else {
485 real_type v (u / scale);
486 sum_squares += v * v;
487 }
488 ++ it;
489 }
490 return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares));
491#endif
492 }
493 // Sparse case
494 template<class I>
495 static BOOST_UBLAS_INLINE
496 result_type apply (I it, const I &it_end) {
497#ifndef BOOST_UBLAS_SCALED_NORM
498 real_type t = real_type ();
499 while (it != it_end) {
500 real_type u (type_traits<value_type>::norm_2 (*it));
501 t += u * u;
502 ++ it;
503 }
504 return static_cast<result_type>(type_traits<real_type>::type_sqrt (t));
505#else
506 real_type scale = real_type ();
507 real_type sum_squares (1);
508 while (it != it_end) {
509 real_type u (type_traits<value_type>::norm_2 (*it));
510 if (scale < u) {
511 real_type v (scale / u);
512 sum_squares = sum_squares * v * v + real_type (1);
513 scale = u;
514 } else {
515 real_type v (u / scale);
516 sum_squares += v * v;
517 }
518 ++ it;
519 }
520 return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares));
521#endif
522 }
523 };
524
525 template<class V>
526 struct vector_norm_2_square :
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_2 (e () (i)));
540 t += u * u;
541 }
542 return t;
543 }
544 // Dense case
545 template<class D, class I>
546 static BOOST_UBLAS_INLINE
547 result_type apply (D size, I it) {
548 real_type t = real_type ();
549 while (-- size >= 0) {
550 real_type u (type_traits<value_type>::norm_2 (*it));
551 t += u * u;
552 ++ it;
553 }
554 return t;
555 }
556 // Sparse case
557 template<class I>
558 static BOOST_UBLAS_INLINE
559 result_type apply (I it, const I &it_end) {
560 real_type t = real_type ();
561 while (it != it_end) {
562 real_type u (type_traits<value_type>::norm_2 (*it));
563 t += u * u;
564 ++ it;
565 }
566 return t;
567 }
568 };
569
570 template<class V>
571 struct vector_norm_inf:
572 public vector_scalar_real_unary_functor<V> {
573 typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
574 typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
575 typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
576
577 template<class E>
578 static BOOST_UBLAS_INLINE
579 result_type apply (const vector_expression<E> &e) {
580 real_type t = real_type ();
581 typedef typename E::size_type vector_size_type;
582 vector_size_type size (e ().size ());
583 for (vector_size_type i = 0; i < size; ++ i) {
584 real_type u (type_traits<value_type>::norm_inf (e () (i)));
585 if (u > t)
586 t = u;
587 }
588 return t;
589 }
590 // Dense case
591 template<class D, class I>
592 static BOOST_UBLAS_INLINE
593 result_type apply (D size, I it) {
594 real_type t = real_type ();
595 while (-- size >= 0) {
596 real_type u (type_traits<value_type>::norm_inf (*it));
597 if (u > t)
598 t = u;
599 ++ it;
600 }
601 return t;
602 }
603 // Sparse case
604 template<class I>
605 static BOOST_UBLAS_INLINE
606 result_type apply (I it, const I &it_end) {
607 real_type t = real_type ();
608 while (it != it_end) {
609 real_type u (type_traits<value_type>::norm_inf (*it));
610 if (u > t)
611 t = u;
612 ++ it;
613 }
614 return t;
615 }
616 };
617
618 // Unary returning index
619 template<class V>
620 struct vector_scalar_index_unary_functor {
621 typedef typename V::value_type value_type;
622 typedef typename type_traits<value_type>::real_type real_type;
623 typedef typename V::size_type result_type;
624 };
625
626 template<class V>
627 struct vector_index_norm_inf:
628 public vector_scalar_index_unary_functor<V> {
629 typedef typename vector_scalar_index_unary_functor<V>::value_type value_type;
630 typedef typename vector_scalar_index_unary_functor<V>::real_type real_type;
631 typedef typename vector_scalar_index_unary_functor<V>::result_type result_type;
632
633 template<class E>
634 static BOOST_UBLAS_INLINE
635 result_type apply (const vector_expression<E> &e) {
636 // ISSUE For CBLAS compatibility return 0 index in empty case
637 result_type i_norm_inf (0);
638 real_type t = real_type ();
639 typedef typename E::size_type vector_size_type;
640 vector_size_type size (e ().size ());
641 for (vector_size_type i = 0; i < size; ++ i) {
642 real_type u (type_traits<value_type>::norm_inf (e () (i)));
643 if (u > t) {
644 i_norm_inf = i;
645 t = u;
646 }
647 }
648 return i_norm_inf;
649 }
650 // Dense case
651 template<class D, class I>
652 static BOOST_UBLAS_INLINE
653 result_type apply (D size, I it) {
654 // ISSUE For CBLAS compatibility return 0 index in empty case
655 result_type i_norm_inf (0);
656 real_type t = real_type ();
657 while (-- size >= 0) {
658 real_type u (type_traits<value_type>::norm_inf (*it));
659 if (u > t) {
660 i_norm_inf = it.index ();
661 t = u;
662 }
663 ++ it;
664 }
665 return i_norm_inf;
666 }
667 // Sparse case
668 template<class I>
669 static BOOST_UBLAS_INLINE
670 result_type apply (I it, const I &it_end) {
671 // ISSUE For CBLAS compatibility return 0 index in empty case
672 result_type i_norm_inf (0);
673 real_type t = real_type ();
674 while (it != it_end) {
675 real_type u (type_traits<value_type>::norm_inf (*it));
676 if (u > t) {
677 i_norm_inf = it.index ();
678 t = u;
679 }
680 ++ it;
681 }
682 return i_norm_inf;
683 }
684 };
685
686 // Binary returning scalar
687 template<class V1, class V2, class TV>
688 struct vector_scalar_binary_functor {
689 typedef TV value_type;
690 typedef TV result_type;
691 };
692
693 template<class V1, class V2, class TV>
694 struct vector_inner_prod:
695 public vector_scalar_binary_functor<V1, V2, TV> {
696 typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type;
697 typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type;
698
699 template<class C1, class C2>
700 static BOOST_UBLAS_INLINE
701 result_type apply (const vector_container<C1> &c1,
702 const vector_container<C2> &c2) {
703#ifdef BOOST_UBLAS_USE_SIMD
704 using namespace raw;
705 typedef typename C1::size_type vector_size_type;
706 vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
707 const typename V1::value_type *data1 = data_const (c1 ());
708 const typename V1::value_type *data2 = data_const (c2 ());
709 vector_size_type s1 = stride (c1 ());
710 vector_size_type s2 = stride (c2 ());
711 result_type t = result_type (0);
712 if (s1 == 1 && s2 == 1) {
713 for (vector_size_type i = 0; i < size; ++ i)
714 t += data1 [i] * data2 [i];
715 } else if (s2 == 1) {
716 for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
717 t += data1 [i1] * data2 [i];
718 } else if (s1 == 1) {
719 for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
720 t += data1 [i] * data2 [i2];
721 } else {
722 for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
723 t += data1 [i1] * data2 [i2];
724 }
725 return t;
726#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
727 return boost::numeric::bindings::atlas::dot (c1 (), c2 ());
728#else
729 return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2));
730#endif
731 }
732 template<class E1, class E2>
733 static BOOST_UBLAS_INLINE
734 result_type apply (const vector_expression<E1> &e1,
735 const vector_expression<E2> &e2) {
736 typedef typename E1::size_type vector_size_type;
737 vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
738 result_type t = result_type (0);
739#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
740 for (vector_size_type i = 0; i < size; ++ i)
741 t += e1 () (i) * e2 () (i);
742#else
743 vector_size_type i (0);
744 DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));
745#endif
746 return t;
747 }
748 // Dense case
749 template<class D, class I1, class I2>
750 static BOOST_UBLAS_INLINE
751 result_type apply (D size, I1 it1, I2 it2) {
752 result_type t = result_type (0);
753#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
754 while (-- size >= 0)
755 t += *it1 * *it2, ++ it1, ++ it2;
756#else
757 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
758#endif
759 return t;
760 }
761 // Packed case
762 template<class I1, class I2>
763 static BOOST_UBLAS_INLINE
764 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
765 result_type t = result_type (0);
766 typedef typename I1::difference_type vector_difference_type;
767 vector_difference_type it1_size (it1_end - it1);
768 vector_difference_type it2_size (it2_end - it2);
769 vector_difference_type diff (0);
770 if (it1_size > 0 && it2_size > 0)
771 diff = it2.index () - it1.index ();
772 if (diff != 0) {
773 vector_difference_type size = (std::min) (diff, it1_size);
774 if (size > 0) {
775 it1 += size;
776 it1_size -= size;
777 diff -= size;
778 }
779 size = (std::min) (- diff, it2_size);
780 if (size > 0) {
781 it2 += size;
782 it2_size -= size;
783 diff += size;
784 }
785 }
786 vector_difference_type size ((std::min) (it1_size, it2_size));
787 while (-- size >= 0)
788 t += *it1 * *it2, ++ it1, ++ it2;
789 return t;
790 }
791 // Sparse case
792 template<class I1, class I2>
793 static BOOST_UBLAS_INLINE
794 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
795 result_type t = result_type (0);
796 if (it1 != it1_end && it2 != it2_end) {
797 for (;;) {
798 if (it1.index () == it2.index ()) {
799 t += *it1 * *it2, ++ it1, ++ it2;
800 if (it1 == it1_end || it2 == it2_end)
801 break;
802 } else if (it1.index () < it2.index ()) {
803 increment (it1, it1_end, it2.index () - it1.index ());
804 if (it1 == it1_end)
805 break;
806 } else if (it1.index () > it2.index ()) {
807 increment (it2, it2_end, it1.index () - it2.index ());
808 if (it2 == it2_end)
809 break;
810 }
811 }
812 }
813 return t;
814 }
815 };
816
817 // Matrix functors
818
819 // Binary returning vector
820 template<class M1, class M2, class TV>
821 struct matrix_vector_binary_functor {
822 typedef typename M1::size_type size_type;
823 typedef typename M1::difference_type difference_type;
824 typedef TV value_type;
825 typedef TV result_type;
826 };
827
828 template<class M1, class M2, class TV>
829 struct matrix_vector_prod1:
830 public matrix_vector_binary_functor<M1, M2, TV> {
831 typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
832 typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
833 typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
834 typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;
835
836 template<class C1, class C2>
837 static BOOST_UBLAS_INLINE
838 result_type apply (const matrix_container<C1> &c1,
839 const vector_container<C2> &c2,
840 size_type i) {
841#ifdef BOOST_UBLAS_USE_SIMD
842 using namespace raw;
843 size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ());
844 const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
845 const typename M2::value_type *data2 = data_const (c2 ());
846 size_type s1 = stride2 (c1 ());
847 size_type s2 = stride (c2 ());
848 result_type t = result_type (0);
849 if (s1 == 1 && s2 == 1) {
850 for (size_type j = 0; j < size; ++ j)
851 t += data1 [j] * data2 [j];
852 } else if (s2 == 1) {
853 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
854 t += data1 [j1] * data2 [j];
855 } else if (s1 == 1) {
856 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
857 t += data1 [j] * data2 [j2];
858 } else {
859 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
860 t += data1 [j1] * data2 [j2];
861 }
862 return t;
863#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
864 return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ());
865#else
866 return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i));
867#endif
868 }
869 template<class E1, class E2>
870 static BOOST_UBLAS_INLINE
871 result_type apply (const matrix_expression<E1> &e1,
872 const vector_expression<E2> &e2,
873 size_type i) {
874 size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
875 result_type t = result_type (0);
876#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
877 for (size_type j = 0; j < size; ++ j)
878 t += e1 () (i, j) * e2 () (j);
879#else
880 size_type j (0);
881 DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j));
882#endif
883 return t;
884 }
885 // Dense case
886 template<class I1, class I2>
887 static BOOST_UBLAS_INLINE
888 result_type apply (difference_type size, I1 it1, I2 it2) {
889 result_type t = result_type (0);
890#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
891 while (-- size >= 0)
892 t += *it1 * *it2, ++ it1, ++ it2;
893#else
894 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
895#endif
896 return t;
897 }
898 // Packed case
899 template<class I1, class I2>
900 static BOOST_UBLAS_INLINE
901 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
902 result_type t = result_type (0);
903 difference_type it1_size (it1_end - it1);
904 difference_type it2_size (it2_end - it2);
905 difference_type diff (0);
906 if (it1_size > 0 && it2_size > 0)
907 diff = it2.index () - it1.index2 ();
908 if (diff != 0) {
909 difference_type size = (std::min) (diff, it1_size);
910 if (size > 0) {
911 it1 += size;
912 it1_size -= size;
913 diff -= size;
914 }
915 size = (std::min) (- diff, it2_size);
916 if (size > 0) {
917 it2 += size;
918 it2_size -= size;
919 diff += size;
920 }
921 }
922 difference_type size ((std::min) (it1_size, it2_size));
923 while (-- size >= 0)
924 t += *it1 * *it2, ++ it1, ++ it2;
925 return t;
926 }
927 // Sparse case
928 template<class I1, class I2>
929 static BOOST_UBLAS_INLINE
930 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
931 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
932 result_type t = result_type (0);
933 if (it1 != it1_end && it2 != it2_end) {
934 size_type it1_index = it1.index2 (), it2_index = it2.index ();
935 for (;;) {
936 difference_type compare = it1_index - it2_index;
937 if (compare == 0) {
938 t += *it1 * *it2, ++ it1, ++ it2;
939 if (it1 != it1_end && it2 != it2_end) {
940 it1_index = it1.index2 ();
941 it2_index = it2.index ();
942 } else
943 break;
944 } else if (compare < 0) {
945 increment (it1, it1_end, - compare);
946 if (it1 != it1_end)
947 it1_index = it1.index2 ();
948 else
949 break;
950 } else if (compare > 0) {
951 increment (it2, it2_end, compare);
952 if (it2 != it2_end)
953 it2_index = it2.index ();
954 else
955 break;
956 }
957 }
958 }
959 return t;
960 }
961 // Sparse packed case
962 template<class I1, class I2>
963 static BOOST_UBLAS_INLINE
964 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
965 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
966 result_type t = result_type (0);
967 while (it1 != it1_end) {
968 t += *it1 * it2 () (it1.index2 ());
969 ++ it1;
970 }
971 return t;
972 }
973 // Packed sparse case
974 template<class I1, class I2>
975 static BOOST_UBLAS_INLINE
976 result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
977 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
978 result_type t = result_type (0);
979 while (it2 != it2_end) {
980 t += it1 () (it1.index1 (), it2.index ()) * *it2;
981 ++ it2;
982 }
983 return t;
984 }
985 // Another dispatcher
986 template<class I1, class I2>
987 static BOOST_UBLAS_INLINE
988 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
989 sparse_bidirectional_iterator_tag) {
990 typedef typename I1::iterator_category iterator1_category;
991 typedef typename I2::iterator_category iterator2_category;
992 return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
993 }
994 };
995
996 template<class M1, class M2, class TV>
997 struct matrix_vector_prod2:
998 public matrix_vector_binary_functor<M1, M2, TV> {
999 typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
1000 typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
1001 typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
1002 typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;
1003
1004 template<class C1, class C2>
1005 static BOOST_UBLAS_INLINE
1006 result_type apply (const vector_container<C1> &c1,
1007 const matrix_container<C2> &c2,
1008 size_type i) {
1009#ifdef BOOST_UBLAS_USE_SIMD
1010 using namespace raw;
1011 size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ());
1012 const typename M1::value_type *data1 = data_const (c1 ());
1013 const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ());
1014 size_type s1 = stride (c1 ());
1015 size_type s2 = stride1 (c2 ());
1016 result_type t = result_type (0);
1017 if (s1 == 1 && s2 == 1) {
1018 for (size_type j = 0; j < size; ++ j)
1019 t += data1 [j] * data2 [j];
1020 } else if (s2 == 1) {
1021 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
1022 t += data1 [j1] * data2 [j];
1023 } else if (s1 == 1) {
1024 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
1025 t += data1 [j] * data2 [j2];
1026 } else {
1027 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
1028 t += data1 [j1] * data2 [j2];
1029 }
1030 return t;
1031#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1032 return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i));
1033#else
1034 return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
1035#endif
1036 }
1037 template<class E1, class E2>
1038 static BOOST_UBLAS_INLINE
1039 result_type apply (const vector_expression<E1> &e1,
1040 const matrix_expression<E2> &e2,
1041 size_type i) {
1042 size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
1043 result_type t = result_type (0);
1044#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1045 for (size_type j = 0; j < size; ++ j)
1046 t += e1 () (j) * e2 () (j, i);
1047#else
1048 size_type j (0);
1049 DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j));
1050#endif
1051 return t;
1052 }
1053 // Dense case
1054 template<class I1, class I2>
1055 static BOOST_UBLAS_INLINE
1056 result_type apply (difference_type size, I1 it1, I2 it2) {
1057 result_type t = result_type (0);
1058#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1059 while (-- size >= 0)
1060 t += *it1 * *it2, ++ it1, ++ it2;
1061#else
1062 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
1063#endif
1064 return t;
1065 }
1066 // Packed case
1067 template<class I1, class I2>
1068 static BOOST_UBLAS_INLINE
1069 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
1070 result_type t = result_type (0);
1071 difference_type it1_size (it1_end - it1);
1072 difference_type it2_size (it2_end - it2);
1073 difference_type diff (0);
1074 if (it1_size > 0 && it2_size > 0)
1075 diff = it2.index1 () - it1.index ();
1076 if (diff != 0) {
1077 difference_type size = (std::min) (diff, it1_size);
1078 if (size > 0) {
1079 it1 += size;
1080 it1_size -= size;
1081 diff -= size;
1082 }
1083 size = (std::min) (- diff, it2_size);
1084 if (size > 0) {
1085 it2 += size;
1086 it2_size -= size;
1087 diff += size;
1088 }
1089 }
1090 difference_type size ((std::min) (it1_size, it2_size));
1091 while (-- size >= 0)
1092 t += *it1 * *it2, ++ it1, ++ it2;
1093 return t;
1094 }
1095 // Sparse case
1096 template<class I1, class I2>
1097 static BOOST_UBLAS_INLINE
1098 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1099 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
1100 result_type t = result_type (0);
1101 if (it1 != it1_end && it2 != it2_end) {
1102 size_type it1_index = it1.index (), it2_index = it2.index1 ();
1103 for (;;) {
1104 difference_type compare = it1_index - it2_index;
1105 if (compare == 0) {
1106 t += *it1 * *it2, ++ it1, ++ it2;
1107 if (it1 != it1_end && it2 != it2_end) {
1108 it1_index = it1.index ();
1109 it2_index = it2.index1 ();
1110 } else
1111 break;
1112 } else if (compare < 0) {
1113 increment (it1, it1_end, - compare);
1114 if (it1 != it1_end)
1115 it1_index = it1.index ();
1116 else
1117 break;
1118 } else if (compare > 0) {
1119 increment (it2, it2_end, compare);
1120 if (it2 != it2_end)
1121 it2_index = it2.index1 ();
1122 else
1123 break;
1124 }
1125 }
1126 }
1127 return t;
1128 }
1129 // Packed sparse case
1130 template<class I1, class I2>
1131 static BOOST_UBLAS_INLINE
1132 result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
1133 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
1134 result_type t = result_type (0);
1135 while (it2 != it2_end) {
1136 t += it1 () (it2.index1 ()) * *it2;
1137 ++ it2;
1138 }
1139 return t;
1140 }
1141 // Sparse packed case
1142 template<class I1, class I2>
1143 static BOOST_UBLAS_INLINE
1144 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
1145 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
1146 result_type t = result_type (0);
1147 while (it1 != it1_end) {
1148 t += *it1 * it2 () (it1.index (), it2.index2 ());
1149 ++ it1;
1150 }
1151 return t;
1152 }
1153 // Another dispatcher
1154 template<class I1, class I2>
1155 static BOOST_UBLAS_INLINE
1156 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1157 sparse_bidirectional_iterator_tag) {
1158 typedef typename I1::iterator_category iterator1_category;
1159 typedef typename I2::iterator_category iterator2_category;
1160 return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
1161 }
1162 };
1163
1164 // Binary returning matrix
1165 template<class M1, class M2, class TV>
1166 struct matrix_matrix_binary_functor {
1167 typedef typename M1::size_type size_type;
1168 typedef typename M1::difference_type difference_type;
1169 typedef TV value_type;
1170 typedef TV result_type;
1171 };
1172
1173 template<class M1, class M2, class TV>
1174 struct matrix_matrix_prod:
1175 public matrix_matrix_binary_functor<M1, M2, TV> {
1176 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::size_type size_type;
1177 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::difference_type difference_type;
1178 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::value_type value_type;
1179 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::result_type result_type;
1180
1181 template<class C1, class C2>
1182 static BOOST_UBLAS_INLINE
1183 result_type apply (const matrix_container<C1> &c1,
1184 const matrix_container<C2> &c2,
1185 size_type i, size_type j) {
1186#ifdef BOOST_UBLAS_USE_SIMD
1187 using namespace raw;
1188 size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ());
1189 const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
1190 const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ());
1191 size_type s1 = stride2 (c1 ());
1192 size_type s2 = stride1 (c2 ());
1193 result_type t = result_type (0);
1194 if (s1 == 1 && s2 == 1) {
1195 for (size_type k = 0; k < size; ++ k)
1196 t += data1 [k] * data2 [k];
1197 } else if (s2 == 1) {
1198 for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1)
1199 t += data1 [k1] * data2 [k];
1200 } else if (s1 == 1) {
1201 for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2)
1202 t += data1 [k] * data2 [k2];
1203 } else {
1204 for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2)
1205 t += data1 [k1] * data2 [k2];
1206 }
1207 return t;
1208#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1209 return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j));
1210#else
1211 boost::ignore_unused(j);
1212 return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
1213#endif
1214 }
1215 template<class E1, class E2>
1216 static BOOST_UBLAS_INLINE
1217 result_type apply (const matrix_expression<E1> &e1,
1218 const matrix_expression<E2> &e2,
1219 size_type i, size_type j) {
1220 size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
1221 result_type t = result_type (0);
1222#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1223 for (size_type k = 0; k < size; ++ k)
1224 t += e1 () (i, k) * e2 () (k, j);
1225#else
1226 size_type k (0);
1227 DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k));
1228#endif
1229 return t;
1230 }
1231 // Dense case
1232 template<class I1, class I2>
1233 static BOOST_UBLAS_INLINE
1234 result_type apply (difference_type size, I1 it1, I2 it2) {
1235 result_type t = result_type (0);
1236#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1237 while (-- size >= 0)
1238 t += *it1 * *it2, ++ it1, ++ it2;
1239#else
1240 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
1241#endif
1242 return t;
1243 }
1244 // Packed case
1245 template<class I1, class I2>
1246 static BOOST_UBLAS_INLINE
1247 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) {
1248 result_type t = result_type (0);
1249 difference_type it1_size (it1_end - it1);
1250 difference_type it2_size (it2_end - it2);
1251 difference_type diff (0);
1252 if (it1_size > 0 && it2_size > 0)
1253 diff = it2.index1 () - it1.index2 ();
1254 if (diff != 0) {
1255 difference_type size = (std::min) (diff, it1_size);
1256 if (size > 0) {
1257 it1 += size;
1258 it1_size -= size;
1259 diff -= size;
1260 }
1261 size = (std::min) (- diff, it2_size);
1262 if (size > 0) {
1263 it2 += size;
1264 it2_size -= size;
1265 diff += size;
1266 }
1267 }
1268 difference_type size ((std::min) (it1_size, it2_size));
1269 while (-- size >= 0)
1270 t += *it1 * *it2, ++ it1, ++ it2;
1271 return t;
1272 }
1273 // Sparse case
1274 template<class I1, class I2>
1275 static BOOST_UBLAS_INLINE
1276 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
1277 result_type t = result_type (0);
1278 if (it1 != it1_end && it2 != it2_end) {
1279 size_type it1_index = it1.index2 (), it2_index = it2.index1 ();
1280 for (;;) {
1281 difference_type compare = difference_type (it1_index - it2_index);
1282 if (compare == 0) {
1283 t += *it1 * *it2, ++ it1, ++ it2;
1284 if (it1 != it1_end && it2 != it2_end) {
1285 it1_index = it1.index2 ();
1286 it2_index = it2.index1 ();
1287 } else
1288 break;
1289 } else if (compare < 0) {
1290 increment (it1, it1_end, - compare);
1291 if (it1 != it1_end)
1292 it1_index = it1.index2 ();
1293 else
1294 break;
1295 } else if (compare > 0) {
1296 increment (it2, it2_end, compare);
1297 if (it2 != it2_end)
1298 it2_index = it2.index1 ();
1299 else
1300 break;
1301 }
1302 }
1303 }
1304 return t;
1305 }
1306 };
1307
1308 // Unary returning scalar norm
1309 template<class M>
1310 struct matrix_scalar_real_unary_functor {
1311 typedef typename M::value_type value_type;
1312 typedef typename type_traits<value_type>::real_type real_type;
1313 typedef real_type result_type;
1314 };
1315
1316 template<class M>
1317 struct matrix_norm_1:
1318 public matrix_scalar_real_unary_functor<M> {
1319 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
1320 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
1321 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
1322
1323 template<class E>
1324 static BOOST_UBLAS_INLINE
1325 result_type apply (const matrix_expression<E> &e) {
1326 real_type t = real_type ();
1327 typedef typename E::size_type matrix_size_type;
1328 matrix_size_type size2 (e ().size2 ());
1329 for (matrix_size_type j = 0; j < size2; ++ j) {
1330 real_type u = real_type ();
1331 matrix_size_type size1 (e ().size1 ());
1332 for (matrix_size_type i = 0; i < size1; ++ i) {
1333 real_type v (type_traits<value_type>::norm_1 (e () (i, j)));
1334 u += v;
1335 }
1336 if (u > t)
1337 t = u;
1338 }
1339 return t;
1340 }
1341 };
1342
1343 template<class M>
1344 struct matrix_norm_frobenius:
1345 public matrix_scalar_real_unary_functor<M> {
1346 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
1347 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
1348 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
1349
1350 template<class E>
1351 static BOOST_UBLAS_INLINE
1352 result_type apply (const matrix_expression<E> &e) {
1353 real_type t = real_type ();
1354 typedef typename E::size_type matrix_size_type;
1355 matrix_size_type size1 (e ().size1 ());
1356 for (matrix_size_type i = 0; i < size1; ++ i) {
1357 matrix_size_type size2 (e ().size2 ());
1358 for (matrix_size_type j = 0; j < size2; ++ j) {
1359 real_type u (type_traits<value_type>::norm_2 (e () (i, j)));
1360 t += u * u;
1361 }
1362 }
1363 return type_traits<real_type>::type_sqrt (t);
1364 }
1365 };
1366
1367 template<class M>
1368 struct matrix_norm_inf:
1369 public matrix_scalar_real_unary_functor<M> {
1370 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
1371 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
1372 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
1373
1374 template<class E>
1375 static BOOST_UBLAS_INLINE
1376 result_type apply (const matrix_expression<E> &e) {
1377 real_type t = real_type ();
1378 typedef typename E::size_type matrix_size_type;
1379 matrix_size_type size1 (e ().size1 ());
1380 for (matrix_size_type i = 0; i < size1; ++ i) {
1381 real_type u = real_type ();
1382 matrix_size_type size2 (e ().size2 ());
1383 for (matrix_size_type j = 0; j < size2; ++ j) {
1384 real_type v (type_traits<value_type>::norm_inf (e () (i, j)));
1385 u += v;
1386 }
1387 if (u > t)
1388 t = u;
1389 }
1390 return t;
1391 }
1392 };
1393
1394 // forward declaration
1395 template <class Z, class D> struct basic_column_major;
1396
1397 // This functor defines storage layout and it's properties
1398 // matrix (i,j) -> storage [i * size_i + j]
1399 template <class Z, class D>
1400 struct basic_row_major {
1401 typedef Z size_type;
1402 typedef D difference_type;
1403 typedef row_major_tag orientation_category;
1404 typedef basic_column_major<Z,D> transposed_layout;
1405
1406 static
1407 BOOST_UBLAS_INLINE
1408 size_type storage_size (size_type size_i, size_type size_j) {
1409 // Guard against size_type overflow
1410 BOOST_UBLAS_CHECK (size_j == 0 || size_i <= (std::numeric_limits<size_type>::max) () / size_j, bad_size ());
1411 return size_i * size_j;
1412 }
1413
1414 // Indexing conversion to storage element
1415 static
1416 BOOST_UBLAS_INLINE
1417 size_type element (size_type i, size_type size_i, size_type j, size_type size_j) {
1418 BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1419 BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1420 detail::ignore_unused_variable_warning(size_i);
1421 // Guard against size_type overflow
1422 BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ());
1423 return i * size_j + j;
1424 }
1425 static
1426 BOOST_UBLAS_INLINE
1427 size_type address (size_type i, size_type size_i, size_type j, size_type size_j) {
1428 BOOST_UBLAS_CHECK (i <= size_i, bad_index ());
1429 BOOST_UBLAS_CHECK (j <= size_j, bad_index ());
1430 // Guard against size_type overflow - address may be size_j past end of storage
1431 BOOST_UBLAS_CHECK (size_j == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ());
1432 detail::ignore_unused_variable_warning(size_i);
1433 return i * size_j + j;
1434 }
1435
1436 // Storage element to index conversion
1437 static
1438 BOOST_UBLAS_INLINE
1439 difference_type distance_i (difference_type k, size_type /* size_i */, size_type size_j) {
1440 return size_j != 0 ? k / size_j : 0;
1441 }
1442 static
1443 BOOST_UBLAS_INLINE
1444 difference_type distance_j (difference_type k, size_type /* size_i */, size_type /* size_j */) {
1445 return k;
1446 }
1447 static
1448 BOOST_UBLAS_INLINE
1449 size_type index_i (difference_type k, size_type /* size_i */, size_type size_j) {
1450 return size_j != 0 ? k / size_j : 0;
1451 }
1452 static
1453 BOOST_UBLAS_INLINE
1454 size_type index_j (difference_type k, size_type /* size_i */, size_type size_j) {
1455 return size_j != 0 ? k % size_j : 0;
1456 }
1457 static
1458 BOOST_UBLAS_INLINE
1459 bool fast_i () {
1460 return false;
1461 }
1462 static
1463 BOOST_UBLAS_INLINE
1464 bool fast_j () {
1465 return true;
1466 }
1467
1468 // Iterating storage elements
1469 template<class I>
1470 static
1471 BOOST_UBLAS_INLINE
1472 void increment_i (I &it, size_type /* size_i */, size_type size_j) {
1473 it += size_j;
1474 }
1475 template<class I>
1476 static
1477 BOOST_UBLAS_INLINE
1478 void increment_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) {
1479 it += n * size_j;
1480 }
1481 template<class I>
1482 static
1483 BOOST_UBLAS_INLINE
1484 void decrement_i (I &it, size_type /* size_i */, size_type size_j) {
1485 it -= size_j;
1486 }
1487 template<class I>
1488 static
1489 BOOST_UBLAS_INLINE
1490 void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) {
1491 it -= n * size_j;
1492 }
1493 template<class I>
1494 static
1495 BOOST_UBLAS_INLINE
1496 void increment_j (I &it, size_type /* size_i */, size_type /* size_j */) {
1497 ++ it;
1498 }
1499 template<class I>
1500 static
1501 BOOST_UBLAS_INLINE
1502 void increment_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
1503 it += n;
1504 }
1505 template<class I>
1506 static
1507 BOOST_UBLAS_INLINE
1508 void decrement_j (I &it, size_type /* size_i */, size_type /* size_j */) {
1509 -- it;
1510 }
1511 template<class I>
1512 static
1513 BOOST_UBLAS_INLINE
1514 void decrement_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
1515 it -= n;
1516 }
1517
1518 // Triangular access
1519 static
1520 BOOST_UBLAS_INLINE
1521 size_type triangular_size (size_type size_i, size_type size_j) {
1522 size_type size = (std::max) (size_i, size_j);
1523 // Guard against size_type overflow - simplified
1524 BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
1525 return ((size + 1) * size) / 2;
1526 }
1527 static
1528 BOOST_UBLAS_INLINE
1529 size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) {
1530 BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1531 BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1532 BOOST_UBLAS_CHECK (i >= j, bad_index ());
1533 detail::ignore_unused_variable_warning(size_i);
1534 detail::ignore_unused_variable_warning(size_j);
1535 // FIXME size_type overflow
1536 // sigma_i (i + 1) = (i + 1) * i / 2
1537 // i = 0 1 2 3, sigma = 0 1 3 6
1538 return ((i + 1) * i) / 2 + j;
1539 }
1540 static
1541 BOOST_UBLAS_INLINE
1542 size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) {
1543 BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1544 BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1545 BOOST_UBLAS_CHECK (i <= j, bad_index ());
1546 // FIXME size_type overflow
1547 // sigma_i (size - i) = size * i - i * (i - 1) / 2
1548 // i = 0 1 2 3, sigma = 0 4 7 9
1549 return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i;
1550 }
1551
1552 // Major and minor indices
1553 static
1554 BOOST_UBLAS_INLINE
1555 size_type index_M (size_type index1, size_type /* index2 */) {
1556 return index1;
1557 }
1558 static
1559 BOOST_UBLAS_INLINE
1560 size_type index_m (size_type /* index1 */, size_type index2) {
1561 return index2;
1562 }
1563 static
1564 BOOST_UBLAS_INLINE
1565 size_type size_M (size_type size_i, size_type /* size_j */) {
1566 return size_i;
1567 }
1568 static
1569 BOOST_UBLAS_INLINE
1570 size_type size_m (size_type /* size_i */, size_type size_j) {
1571 return size_j;
1572 }
1573 };
1574
1575 // This functor defines storage layout and it's properties
1576 // matrix (i,j) -> storage [i + j * size_i]
1577 template <class Z, class D>
1578 struct basic_column_major {
1579 typedef Z size_type;
1580 typedef D difference_type;
1581 typedef column_major_tag orientation_category;
1582 typedef basic_row_major<Z,D> transposed_layout;
1583
1584 static
1585 BOOST_UBLAS_INLINE
1586 size_type storage_size (size_type size_i, size_type size_j) {
1587 // Guard against size_type overflow
1588 BOOST_UBLAS_CHECK (size_i == 0 || size_j <= (std::numeric_limits<size_type>::max) () / size_i, bad_size ());
1589 return size_i * size_j;
1590 }
1591
1592 // Indexing conversion to storage element
1593 static
1594 BOOST_UBLAS_INLINE
1595 size_type element (size_type i, size_type size_i, size_type j, size_type size_j) {
1596 BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1597 BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1598 detail::ignore_unused_variable_warning(size_j);
1599 // Guard against size_type overflow
1600 BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ());
1601 return i + j * size_i;
1602 }
1603 static
1604 BOOST_UBLAS_INLINE
1605 size_type address (size_type i, size_type size_i, size_type j, size_type size_j) {
1606 BOOST_UBLAS_CHECK (i <= size_i, bad_index ());
1607 BOOST_UBLAS_CHECK (j <= size_j, bad_index ());
1608 detail::ignore_unused_variable_warning(size_j);
1609 // Guard against size_type overflow - address may be size_i past end of storage
1610 BOOST_UBLAS_CHECK (size_i == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ());
1611 return i + j * size_i;
1612 }
1613
1614 // Storage element to index conversion
1615 static
1616 BOOST_UBLAS_INLINE
1617 difference_type distance_i (difference_type k, size_type /* size_i */, size_type /* size_j */) {
1618 return k;
1619 }
1620 static
1621 BOOST_UBLAS_INLINE
1622 difference_type distance_j (difference_type k, size_type size_i, size_type /* size_j */) {
1623 return size_i != 0 ? k / size_i : 0;
1624 }
1625 static
1626 BOOST_UBLAS_INLINE
1627 size_type index_i (difference_type k, size_type size_i, size_type /* size_j */) {
1628 return size_i != 0 ? k % size_i : 0;
1629 }
1630 static
1631 BOOST_UBLAS_INLINE
1632 size_type index_j (difference_type k, size_type size_i, size_type /* size_j */) {
1633 return size_i != 0 ? k / size_i : 0;
1634 }
1635 static
1636 BOOST_UBLAS_INLINE
1637 bool fast_i () {
1638 return true;
1639 }
1640 static
1641 BOOST_UBLAS_INLINE
1642 bool fast_j () {
1643 return false;
1644 }
1645
1646 // Iterating
1647 template<class I>
1648 static
1649 BOOST_UBLAS_INLINE
1650 void increment_i (I &it, size_type /* size_i */, size_type /* size_j */) {
1651 ++ it;
1652 }
1653 template<class I>
1654 static
1655 BOOST_UBLAS_INLINE
1656 void increment_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
1657 it += n;
1658 }
1659 template<class I>
1660 static
1661 BOOST_UBLAS_INLINE
1662 void decrement_i (I &it, size_type /* size_i */, size_type /* size_j */) {
1663 -- it;
1664 }
1665 template<class I>
1666 static
1667 BOOST_UBLAS_INLINE
1668 void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
1669 it -= n;
1670 }
1671 template<class I>
1672 static
1673 BOOST_UBLAS_INLINE
1674 void increment_j (I &it, size_type size_i, size_type /* size_j */) {
1675 it += size_i;
1676 }
1677 template<class I>
1678 static
1679 BOOST_UBLAS_INLINE
1680 void increment_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) {
1681 it += n * size_i;
1682 }
1683 template<class I>
1684 static
1685 BOOST_UBLAS_INLINE
1686 void decrement_j (I &it, size_type size_i, size_type /* size_j */) {
1687 it -= size_i;
1688 }
1689 template<class I>
1690 static
1691 BOOST_UBLAS_INLINE
1692 void decrement_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) {
1693 it -= n* size_i;
1694 }
1695
1696 // Triangular access
1697 static
1698 BOOST_UBLAS_INLINE
1699 size_type triangular_size (size_type size_i, size_type size_j) {
1700 size_type size = (std::max) (size_i, size_j);
1701 // Guard against size_type overflow - simplified
1702 BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
1703 return ((size + 1) * size) / 2;
1704 }
1705 static
1706 BOOST_UBLAS_INLINE
1707 size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) {
1708 BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1709 BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1710 BOOST_UBLAS_CHECK (i >= j, bad_index ());
1711 // FIXME size_type overflow
1712 // sigma_j (size - j) = size * j - j * (j - 1) / 2
1713 // j = 0 1 2 3, sigma = 0 4 7 9
1714 return i - j + (j * (2 * (std::max) (size_i, size_j) - j + 1)) / 2;
1715 }
1716 static
1717 BOOST_UBLAS_INLINE
1718 size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) {
1719 BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1720 BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1721 BOOST_UBLAS_CHECK (i <= j, bad_index ());
1722 boost::ignore_unused(size_i, size_j);
1723 // FIXME size_type overflow
1724 // sigma_j (j + 1) = (j + 1) * j / 2
1725 // j = 0 1 2 3, sigma = 0 1 3 6
1726 return i + ((j + 1) * j) / 2;
1727 }
1728
1729 // Major and minor indices
1730 static
1731 BOOST_UBLAS_INLINE
1732 size_type index_M (size_type /* index1 */, size_type index2) {
1733 return index2;
1734 }
1735 static
1736 BOOST_UBLAS_INLINE
1737 size_type index_m (size_type index1, size_type /* index2 */) {
1738 return index1;
1739 }
1740 static
1741 BOOST_UBLAS_INLINE
1742 size_type size_M (size_type /* size_i */, size_type size_j) {
1743 return size_j;
1744 }
1745 static
1746 BOOST_UBLAS_INLINE
1747 size_type size_m (size_type size_i, size_type /* size_j */) {
1748 return size_i;
1749 }
1750 };
1751
1752
1753 template <class Z>
1754 struct basic_full {
1755 typedef Z size_type;
1756
1757 template<class L>
1758 static
1759 BOOST_UBLAS_INLINE
1760 size_type packed_size (L, size_type size_i, size_type size_j) {
1761 return L::storage_size (size_i, size_j);
1762 }
1763
1764 static
1765 BOOST_UBLAS_INLINE
1766 bool zero (size_type /* i */, size_type /* j */) {
1767 return false;
1768 }
1769 static
1770 BOOST_UBLAS_INLINE
1771 bool one (size_type /* i */, size_type /* j */) {
1772 return false;
1773 }
1774 static
1775 BOOST_UBLAS_INLINE
1776 bool other (size_type /* i */, size_type /* j */) {
1777 return true;
1778 }
1779 // FIXME: this should not be used at all
1780 static
1781 BOOST_UBLAS_INLINE
1782 size_type restrict1 (size_type i, size_type /* j */) {
1783 return i;
1784 }
1785 static
1786 BOOST_UBLAS_INLINE
1787 size_type restrict2 (size_type /* i */, size_type j) {
1788 return j;
1789 }
1790 static
1791 BOOST_UBLAS_INLINE
1792 size_type mutable_restrict1 (size_type i, size_type /* j */) {
1793 return i;
1794 }
1795 static
1796 BOOST_UBLAS_INLINE
1797 size_type mutable_restrict2 (size_type /* i */, size_type j) {
1798 return j;
1799 }
1800 };
1801
1802 namespace detail {
1803 template < class L >
1804 struct transposed_structure {
1805 typedef typename L::size_type size_type;
1806
1807 template<class LAYOUT>
1808 static
1809 BOOST_UBLAS_INLINE
1810 size_type packed_size (LAYOUT l, size_type size_i, size_type size_j) {
1811 return L::packed_size(l, size_j, size_i);
1812 }
1813
1814 static
1815 BOOST_UBLAS_INLINE
1816 bool zero (size_type i, size_type j) {
1817 return L::zero(j, i);
1818 }
1819 static
1820 BOOST_UBLAS_INLINE
1821 bool one (size_type i, size_type j) {
1822 return L::one(j, i);
1823 }
1824 static
1825 BOOST_UBLAS_INLINE
1826 bool other (size_type i, size_type j) {
1827 return L::other(j, i);
1828 }
1829 template<class LAYOUT>
1830 static
1831 BOOST_UBLAS_INLINE
1832 size_type element (LAYOUT /* l */, size_type i, size_type size_i, size_type j, size_type size_j) {
1833 return L::element(typename LAYOUT::transposed_layout(), j, size_j, i, size_i);
1834 }
1835
1836 static
1837 BOOST_UBLAS_INLINE
1838 size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
1839 return L::restrict2(j, i, size2, size1);
1840 }
1841 static
1842 BOOST_UBLAS_INLINE
1843 size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
1844 return L::restrict1(j, i, size2, size1);
1845 }
1846 static
1847 BOOST_UBLAS_INLINE
1848 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
1849 return L::mutable_restrict2(j, i, size2, size1);
1850 }
1851 static
1852 BOOST_UBLAS_INLINE
1853 size_type mutable_restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
1854 return L::mutable_restrict1(j, i, size2, size1);
1855 }
1856
1857 static
1858 BOOST_UBLAS_INLINE
1859 size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
1860 return L::global_restrict2(index2, size2, index1, size1);
1861 }
1862 static
1863 BOOST_UBLAS_INLINE
1864 size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
1865 return L::global_restrict1(index2, size2, index1, size1);
1866 }
1867 static
1868 BOOST_UBLAS_INLINE
1869 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
1870 return L::global_mutable_restrict2(index2, size2, index1, size1);
1871 }
1872 static
1873 BOOST_UBLAS_INLINE
1874 size_type global_mutable_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
1875 return L::global_mutable_restrict1(index2, size2, index1, size1);
1876 }
1877 };
1878 }
1879
1880 template <class Z>
1881 struct basic_lower {
1882 typedef Z size_type;
1883 typedef lower_tag triangular_type;
1884
1885 template<class L>
1886 static
1887 BOOST_UBLAS_INLINE
1888 size_type packed_size (L, size_type size_i, size_type size_j) {
1889 return L::triangular_size (size_i, size_j);
1890 }
1891
1892 static
1893 BOOST_UBLAS_INLINE
1894 bool zero (size_type i, size_type j) {
1895 return j > i;
1896 }
1897 static
1898 BOOST_UBLAS_INLINE
1899 bool one (size_type /* i */, size_type /* j */) {
1900 return false;
1901 }
1902 static
1903 BOOST_UBLAS_INLINE
1904 bool other (size_type i, size_type j) {
1905 return j <= i;
1906 }
1907 template<class L>
1908 static
1909 BOOST_UBLAS_INLINE
1910 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
1911 return L::lower_element (i, size_i, j, size_j);
1912 }
1913
1914 // return nearest valid index in column j
1915 static
1916 BOOST_UBLAS_INLINE
1917 size_type restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
1918 return (std::max)(j, (std::min) (size1, i));
1919 }
1920 // return nearest valid index in row i
1921 static
1922 BOOST_UBLAS_INLINE
1923 size_type restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
1924 return (std::max)(size_type(0), (std::min) (i+1, j));
1925 }
1926 // return nearest valid mutable index in column j
1927 static
1928 BOOST_UBLAS_INLINE
1929 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
1930 return (std::max)(j, (std::min) (size1, i));
1931 }
1932 // return nearest valid mutable index in row i
1933 static
1934 BOOST_UBLAS_INLINE
1935 size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
1936 return (std::max)(size_type(0), (std::min) (i+1, j));
1937 }
1938
1939 // return an index between the first and (1+last) filled row
1940 static
1941 BOOST_UBLAS_INLINE
1942 size_type global_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
1943 return (std::max)(size_type(0), (std::min)(size1, index1) );
1944 }
1945 // return an index between the first and (1+last) filled column
1946 static
1947 BOOST_UBLAS_INLINE
1948 size_type global_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
1949 return (std::max)(size_type(0), (std::min)(size2, index2) );
1950 }
1951
1952 // return an index between the first and (1+last) filled mutable row
1953 static
1954 BOOST_UBLAS_INLINE
1955 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
1956 return (std::max)(size_type(0), (std::min)(size1, index1) );
1957 }
1958 // return an index between the first and (1+last) filled mutable column
1959 static
1960 BOOST_UBLAS_INLINE
1961 size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
1962 return (std::max)(size_type(0), (std::min)(size2, index2) );
1963 }
1964 };
1965
1966 // the first row only contains a single 1. Thus it is not stored.
1967 template <class Z>
1968 struct basic_unit_lower : public basic_lower<Z> {
1969 typedef Z size_type;
1970 typedef unit_lower_tag triangular_type;
1971
1972 template<class L>
1973 static
1974 BOOST_UBLAS_INLINE
1975 size_type packed_size (L, size_type size_i, size_type size_j) {
1976 // Zero size strict triangles are bad at this point
1977 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ());
1978 return L::triangular_size (size_i - 1, size_j - 1);
1979 }
1980
1981 static
1982 BOOST_UBLAS_INLINE
1983 bool one (size_type i, size_type j) {
1984 return j == i;
1985 }
1986 static
1987 BOOST_UBLAS_INLINE
1988 bool other (size_type i, size_type j) {
1989 return j < i;
1990 }
1991 template<class L>
1992 static
1993 BOOST_UBLAS_INLINE
1994 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
1995 // Zero size strict triangles are bad at this point
1996 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ());
1997 return L::lower_element (i-1, size_i - 1, j, size_j - 1);
1998 }
1999
2000 static
2001 BOOST_UBLAS_INLINE
2002 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
2003 return (std::max)(j+1, (std::min) (size1, i));
2004 }
2005 static
2006 BOOST_UBLAS_INLINE
2007 size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
2008 return (std::max)(size_type(0), (std::min) (i, j));
2009 }
2010
2011 // return an index between the first and (1+last) filled mutable row
2012 static
2013 BOOST_UBLAS_INLINE
2014 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
2015 return (std::max)(size_type(1), (std::min)(size1, index1) );
2016 }
2017 // return an index between the first and (1+last) filled mutable column
2018 static
2019 BOOST_UBLAS_INLINE
2020 size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
2021 BOOST_UBLAS_CHECK( size2 >= 1 , external_logic() );
2022 return (std::max)(size_type(0), (std::min)(size2-1, index2) );
2023 }
2024 };
2025
2026 // the first row only contains no element. Thus it is not stored.
2027 template <class Z>
2028 struct basic_strict_lower : public basic_unit_lower<Z> {
2029 typedef Z size_type;
2030 typedef strict_lower_tag triangular_type;
2031
2032 template<class L>
2033 static
2034 BOOST_UBLAS_INLINE
2035 size_type packed_size (L, size_type size_i, size_type size_j) {
2036 // Zero size strict triangles are bad at this point
2037 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ());
2038 return L::triangular_size (size_i - 1, size_j - 1);
2039 }
2040
2041 static
2042 BOOST_UBLAS_INLINE
2043 bool zero (size_type i, size_type j) {
2044 return j >= i;
2045 }
2046 static
2047 BOOST_UBLAS_INLINE
2048 bool one (size_type /*i*/, size_type /*j*/) {
2049 return false;
2050 }
2051 static
2052 BOOST_UBLAS_INLINE
2053 bool other (size_type i, size_type j) {
2054 return j < i;
2055 }
2056 template<class L>
2057 static
2058 BOOST_UBLAS_INLINE
2059 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
2060 // Zero size strict triangles are bad at this point
2061 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ());
2062 return L::lower_element (i-1, size_i - 1, j, size_j - 1);
2063 }
2064
2065 static
2066 BOOST_UBLAS_INLINE
2067 size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
2068 return basic_unit_lower<Z>::mutable_restrict1(i, j, size1, size2);
2069 }
2070 static
2071 BOOST_UBLAS_INLINE
2072 size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
2073 return basic_unit_lower<Z>::mutable_restrict2(i, j, size1, size2);
2074 }
2075
2076 // return an index between the first and (1+last) filled row
2077 static
2078 BOOST_UBLAS_INLINE
2079 size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
2080 return basic_unit_lower<Z>::global_mutable_restrict1(index1, size1, index2, size2);
2081 }
2082 // return an index between the first and (1+last) filled column
2083 static
2084 BOOST_UBLAS_INLINE
2085 size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
2086 return basic_unit_lower<Z>::global_mutable_restrict2(index1, size1, index2, size2);
2087 }
2088 };
2089
2090
2091 template <class Z>
2092 struct basic_upper : public detail::transposed_structure<basic_lower<Z> >
2093 {
2094 typedef upper_tag triangular_type;
2095 };
2096
2097 template <class Z>
2098 struct basic_unit_upper : public detail::transposed_structure<basic_unit_lower<Z> >
2099 {
2100 typedef unit_upper_tag triangular_type;
2101 };
2102
2103 template <class Z>
2104 struct basic_strict_upper : public detail::transposed_structure<basic_strict_lower<Z> >
2105 {
2106 typedef strict_upper_tag triangular_type;
2107 };
2108
2109
2110}}}
2111
2112#endif
2113

source code of include/boost/numeric/ublas/functional.hpp