source: NonGTP/Boost/boost/numeric/ublas/functional.hpp @ 857

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