latest from svn
[boost:numeric.git] / include / boost / numeric / ublas / functional.hpp
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
25 namespace 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
39 namespace 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