New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
functional.hpp in vendors/XIOS/current/extern/boost/include/boost/numeric/ublas – NEMO

source: vendors/XIOS/current/extern/boost/include/boost/numeric/ublas/functional.hpp @ 3408

Last change on this file since 3408 was 3408, checked in by rblod, 12 years ago

importing initial XIOS vendor drop

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