- Main Page
- Related Pages
- Classes
- Files

- File List

# functional.hpp

    00001 //
    00002 //  Copyright (c) 2000-2009
    00003 //  Joerg Walter, Mathias Koch, Gunter Winkler
    00004 //
    00005 //  Distributed under the Boost Software License, Version 1.0. (See
    00006 //  accompanying file LICENSE_1_0.txt or copy at
    00007 //  http://www.boost.org/LICENSE_1_0.txt)
    00008 //
    00009 //  The authors gratefully acknowledge the support of
    00010 //  GeNeSys mbH & Co. KG in producing this work.
    00011 //
    00012
    00013 #ifndef _BOOST_UBLAS_FUNCTIONAL_
    00014 #define _BOOST_UBLAS_FUNCTIONAL_
    00015
    00016 #include <functional>
    00017
    00018 #include <boost/numeric/ublas/traits.hpp>
    00019 #ifdef BOOST_UBLAS_USE_DUFF_DEVICE
    00020 #include <boost/numeric/ublas/detail/duff.hpp>
    00021 #endif
    00022 #ifdef BOOST_UBLAS_USE_SIMD
    00023 #include <boost/numeric/ublas/detail/raw.hpp>
    00024 #else
    00025 namespace boost { namespace numeric { namespace ublas { namespace raw {
    00026 }}}}
    00027 #endif
    00028 #ifdef BOOST_UBLAS_HAVE_BINDINGS
    00029 #include <boost/numeric/bindings/traits/std_vector.hpp>
    00030 #include <boost/numeric/bindings/traits/ublas_vector.hpp>
    00031 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
    00032 #include <boost/numeric/bindings/atlas/cblas.hpp>
    00033 #endif
    00034
    00035 #include <boost/numeric/ublas/detail/definitions.hpp>
    00036
    00037
    00038
    00039 namespace boost { namespace numeric { namespace ublas {
    00040
    00041     // Scalar functors
    00042
    00043     // Unary
    00044     template<class T>
    00045     struct scalar_unary_functor {
    00046         typedef T value_type;
    00047         typedef typename type_traits<T>::const_reference argument_type;
    00048         typedef typename type_traits<T>::value_type result_type;
    00049     };
    00050
    00051     template<class T>
    00052     struct scalar_identity:
    00053         public scalar_unary_functor<T> {
    00054         typedef typename scalar_unary_functor<T>::argument_type argument_type;
    00055         typedef typename scalar_unary_functor<T>::result_type result_type;
    00056
    00057         static BOOST_UBLAS_INLINE
    00058         result_type apply (argument_type t) {
    00059             return t;
    00060         }
    00061     };
    00062     template<class T>
    00063     struct scalar_negate:
    00064         public scalar_unary_functor<T> {
    00065         typedef typename scalar_unary_functor<T>::argument_type argument_type;
    00066         typedef typename scalar_unary_functor<T>::result_type result_type;
    00067
    00068         static BOOST_UBLAS_INLINE
    00069         result_type apply (argument_type t) {
    00070             return - t;
    00071         }
    00072     };
    00073     template<class T>
    00074     struct scalar_conj:
    00075         public scalar_unary_functor<T> {
    00076         typedef typename scalar_unary_functor<T>::value_type value_type;
    00077         typedef typename scalar_unary_functor<T>::argument_type argument_type;
    00078         typedef typename scalar_unary_functor<T>::result_type result_type;
    00079
    00080         static BOOST_UBLAS_INLINE
    00081         result_type apply (argument_type t) {
    00082             return type_traits<value_type>::conj (t);
    00083         }
    00084     };
    00085
    00086     // Unary returning real
    00087     template<class T>
    00088     struct scalar_real_unary_functor {
    00089         typedef T value_type;
    00090         typedef typename type_traits<T>::const_reference argument_type;
    00091         typedef typename type_traits<T>::real_type result_type;
    00092     };
    00093
    00094     template<class T>
    00095     struct scalar_real:
    00096         public scalar_real_unary_functor<T> {
    00097         typedef typename scalar_real_unary_functor<T>::value_type value_type;
    00098         typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
    00099         typedef typename scalar_real_unary_functor<T>::result_type result_type;
    00100
    00101         static BOOST_UBLAS_INLINE
    00102         result_type apply (argument_type t) {
    00103             return type_traits<value_type>::real (t);
    00104         }
    00105     };
    00106     template<class T>
    00107     struct scalar_imag:
    00108         public scalar_real_unary_functor<T> {
    00109         typedef typename scalar_real_unary_functor<T>::value_type value_type;
    00110         typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
    00111         typedef typename scalar_real_unary_functor<T>::result_type result_type;
    00112
    00113         static BOOST_UBLAS_INLINE
    00114         result_type apply (argument_type t) {
    00115             return type_traits<value_type>::imag (t);
    00116         }
    00117     };
    00118
    00119     // Binary
    00120     template<class T1, class T2>
    00121     struct scalar_binary_functor {
    00122         typedef typename type_traits<T1>::const_reference argument1_type;
    00123         typedef typename type_traits<T2>::const_reference argument2_type;
    00124         typedef typename promote_traits<T1, T2>::promote_type result_type;
    00125     };
    00126
    00127     template<class T1, class T2>
    00128     struct scalar_plus:
    00129         public scalar_binary_functor<T1, T2> {
    00130         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
    00131         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
    00132         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
    00133
    00134         static BOOST_UBLAS_INLINE
    00135         result_type apply (argument1_type t1, argument2_type t2) {
    00136             return t1 + t2;
    00137         }
    00138     };
    00139     template<class T1, class T2>
    00140     struct scalar_minus:
    00141         public scalar_binary_functor<T1, T2> {
    00142         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
    00143         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
    00144         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
    00145
    00146         static BOOST_UBLAS_INLINE
    00147         result_type apply (argument1_type t1, argument2_type t2) {
    00148             return t1 - t2;
    00149         }
    00150     };
    00151     template<class T1, class T2>
    00152     struct scalar_multiplies:
    00153         public scalar_binary_functor<T1, T2> {
    00154         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
    00155         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
    00156         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
    00157
    00158         static BOOST_UBLAS_INLINE
    00159         result_type apply (argument1_type t1, argument2_type t2) {
    00160             return t1 * t2;
    00161         }
    00162     };
    00163     template<class T1, class T2>
    00164     struct scalar_divides:
    00165         public scalar_binary_functor<T1, T2> {
    00166         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
    00167         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
    00168         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
    00169
    00170         static BOOST_UBLAS_INLINE
    00171         result_type apply (argument1_type t1, argument2_type t2) {
    00172             return t1 / t2;
    00173         }
    00174     };
    00175
    00176     template<class T1, class T2>
    00177     struct scalar_binary_assign_functor {
    00178         // ISSUE Remove reference to avoid reference to reference problems
    00179         typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
    00180         typedef typename type_traits<T2>::const_reference argument2_type;
    00181     };
    00182
    00183     struct assign_tag {};
    00184     struct computed_assign_tag {};
    00185
    00186     template<class T1, class T2>
    00187     struct scalar_assign:
    00188         public scalar_binary_assign_functor<T1, T2> {
    00189         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
    00190         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
    00191 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    00192         static const bool computed ;
    00193 #else
    00194         static const bool computed = false ;
    00195 #endif
    00196
    00197         static BOOST_UBLAS_INLINE
    00198         void apply (argument1_type t1, argument2_type t2) {
    00199             t1 = t2;
    00200         }
    00201
    00202         template<class U1, class U2>
    00203         struct rebind {
    00204             typedef scalar_assign<U1, U2> other;
    00205         };
    00206     };
    00207
    00208 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    00209     template<class T1, class T2>
    00210     const bool scalar_assign<T1,T2>::computed = false;
    00211 #endif
    00212
    00213     template<class T1, class T2>
    00214     struct scalar_plus_assign:
    00215         public scalar_binary_assign_functor<T1, T2> {
    00216         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
    00217         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
    00218 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    00219         static const bool computed ;
    00220 #else
    00221         static const bool computed = true ;
    00222 #endif
    00223
    00224         static BOOST_UBLAS_INLINE
    00225         void apply (argument1_type t1, argument2_type t2) {
    00226             t1 += t2;
    00227         }
    00228
    00229         template<class U1, class U2>
    00230         struct rebind {
    00231             typedef scalar_plus_assign<U1, U2> other;
    00232         };
    00233     };
    00234
    00235 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    00236     template<class T1, class T2>
    00237     const bool scalar_plus_assign<T1,T2>::computed = true;
    00238 #endif
    00239
    00240     template<class T1, class T2>
    00241     struct scalar_minus_assign:
    00242         public scalar_binary_assign_functor<T1, T2> {
    00243         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
    00244         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
    00245 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    00246         static const bool computed ;
    00247 #else
    00248         static const bool computed = true ;
    00249 #endif
    00250
    00251         static BOOST_UBLAS_INLINE
    00252         void apply (argument1_type t1, argument2_type t2) {
    00253             t1 -= t2;
    00254         }
    00255
    00256         template<class U1, class U2>
    00257         struct rebind {
    00258             typedef scalar_minus_assign<U1, U2> other;
    00259         };
    00260     };
    00261
    00262 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    00263     template<class T1, class T2>
    00264     const bool scalar_minus_assign<T1,T2>::computed = true;
    00265 #endif
    00266
    00267     template<class T1, class T2>
    00268     struct scalar_multiplies_assign:
    00269         public scalar_binary_assign_functor<T1, T2> {
    00270         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
    00271         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
    00272         static const bool computed = true;
    00273
    00274         static BOOST_UBLAS_INLINE
    00275         void apply (argument1_type t1, argument2_type t2) {
    00276             t1 *= t2;
    00277         }
    00278
    00279         template<class U1, class U2>
    00280         struct rebind {
    00281             typedef scalar_multiplies_assign<U1, U2> other;
    00282         };
    00283     };
    00284     template<class T1, class T2>
    00285     struct scalar_divides_assign:
    00286         public scalar_binary_assign_functor<T1, T2> {
    00287         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
    00288         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
    00289         static const bool computed ;
    00290
    00291         static BOOST_UBLAS_INLINE
    00292         void apply (argument1_type t1, argument2_type t2) {
    00293             t1 /= t2;
    00294         }
    00295
    00296         template<class U1, class U2>
    00297         struct rebind {
    00298             typedef scalar_divides_assign<U1, U2> other;
    00299         };
    00300     };
    00301     template<class T1, class T2>
    00302     const bool scalar_divides_assign<T1,T2>::computed = true;
    00303
    00304     template<class T1, class T2>
    00305     struct scalar_binary_swap_functor {
    00306         typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
    00307         typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type;
    00308     };
    00309
    00310     template<class T1, class T2>
    00311     struct scalar_swap:
    00312         public scalar_binary_swap_functor<T1, T2> {
    00313         typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type;
    00314         typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type;
    00315
    00316         static BOOST_UBLAS_INLINE
    00317         void apply (argument1_type t1, argument2_type t2) {
    00318             std::swap (t1, t2);
    00319         }
    00320
    00321         template<class U1, class U2>
    00322         struct rebind {
    00323             typedef scalar_swap<U1, U2> other;
    00324         };
    00325     };
    00326
    00327     // Vector functors
    00328
    00329     // Unary returning scalar
    00330     template<class V>
    00331     struct vector_scalar_unary_functor {
    00332         typedef typename V::value_type value_type;
    00333         typedef typename V::value_type result_type;
    00334     };
    00335
    00336     template<class V>
    00337     struct vector_sum:
    00338         public vector_scalar_unary_functor<V> {
    00339         typedef typename vector_scalar_unary_functor<V>::value_type value_type;
    00340         typedef typename vector_scalar_unary_functor<V>::result_type result_type;
    00341
    00342         template<class E>
    00343         static BOOST_UBLAS_INLINE
    00344         result_type apply (const vector_expression<E> &e) {
    00345             result_type t = result_type (0);
    00346             typedef typename E::size_type vector_size_type;
    00347             vector_size_type size (e ().size ());
    00348             for (vector_size_type i = 0; i < size; ++ i)
    00349                 t += e () (i);
    00350             return t;
    00351         }
    00352         // Dense case
    00353         template<class D, class I>
    00354         static BOOST_UBLAS_INLINE
    00355         result_type apply (D size, I it) {
    00356             result_type t = result_type (0);
    00357             while (-- size >= 0)
    00358                 t += *it, ++ it;
    00359             return t;
    00360         }
    00361         // Sparse case
    00362         template<class I>
    00363         static BOOST_UBLAS_INLINE
    00364         result_type apply (I it, const I &it_end) {
    00365             result_type t = result_type (0);
    00366             while (it != it_end)
    00367                 t += *it, ++ it;
    00368             return t;
    00369         }
    00370     };
    00371
    00372     // Unary returning real scalar
    00373     template<class V>
    00374     struct vector_scalar_real_unary_functor {
    00375         typedef typename V::value_type value_type;
    00376         typedef typename type_traits<value_type>::real_type real_type;
    00377         typedef real_type result_type;
    00378     };
    00379
    00380     template<class V>
    00381     struct vector_norm_1:
    00382         public vector_scalar_real_unary_functor<V> {
    00383         typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
    00384         typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
    00385         typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
    00386
    00387         template<class E>
    00388         static BOOST_UBLAS_INLINE
    00389         result_type apply (const vector_expression<E> &e) {
    00390             real_type t = real_type ();
    00391             typedef typename E::size_type vector_size_type;
    00392             vector_size_type size (e ().size ());
    00393             for (vector_size_type i = 0; i < size; ++ i) {
    00394                 real_type u (type_traits<value_type>::type_abs (e () (i)));
    00395                 t += u;
    00396             }
    00397             return t;
    00398         }
    00399         // Dense case
    00400         template<class D, class I>
    00401         static BOOST_UBLAS_INLINE
    00402         result_type apply (D size, I it) {
    00403             real_type t = real_type ();
    00404             while (-- size >= 0) {
    00405                 real_type u (type_traits<value_type>::norm_1 (*it));
    00406                 t += u;
    00407                 ++ it;
    00408             }
    00409             return t;
    00410         }
    00411         // Sparse case
    00412         template<class I>
    00413         static BOOST_UBLAS_INLINE
    00414         result_type apply (I it, const I &it_end) {
    00415             real_type t = real_type ();
    00416             while (it != it_end) {
    00417                 real_type u (type_traits<value_type>::norm_1 (*it));
    00418                 t += u;
    00419                 ++ it;
    00420             }
    00421             return t;
    00422         }
    00423     };
    00424     template<class V>
    00425     struct vector_norm_2:
    00426         public vector_scalar_real_unary_functor<V> {
    00427         typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
    00428         typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
    00429         typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
    00430
    00431         template<class E>
    00432         static BOOST_UBLAS_INLINE
    00433         result_type apply (const vector_expression<E> &e) {
    00434 #ifndef BOOST_UBLAS_SCALED_NORM
    00435             real_type t = real_type ();
    00436             typedef typename E::size_type vector_size_type;
    00437             vector_size_type size (e ().size ());
    00438             for (vector_size_type i = 0; i < size; ++ i) {
    00439                 real_type u (type_traits<value_type>::norm_2 (e () (i)));
    00440                 t +=  u * u;
    00441             }
    00442             return type_traits<real_type>::type_sqrt (t);
    00443 #else
    00444             real_type scale = real_type ();
    00445             real_type sum_squares (1);
    00446             size_type size (e ().size ());
    00447             for (size_type i = 0; i < size; ++ i) {
    00448                 real_type u (type_traits<value_type>::norm_2 (e () (i)));
    00449                 if ( real_type () /* zero */ == u ) continue;
    00450                 if (scale < u) {
    00451                     real_type v (scale / u);
    00452                     sum_squares = sum_squares * v * v + real_type (1);
    00453                     scale = u;
    00454                 } else {
    00455                     real_type v (u / scale);
    00456                     sum_squares += v * v;
    00457                 }
    00458             }
    00459             return scale * type_traits<real_type>::type_sqrt (sum_squares);
    00460 #endif
    00461         }
    00462         // Dense case
    00463         template<class D, class I>
    00464         static BOOST_UBLAS_INLINE
    00465         result_type apply (D size, I it) {
    00466 #ifndef BOOST_UBLAS_SCALED_NORM
    00467             real_type t = real_type ();
    00468             while (-- size >= 0) {
    00469                 real_type u (type_traits<value_type>::norm_2 (*it));
    00470                 t +=  u * u;
    00471                 ++ it;
    00472             }
    00473             return type_traits<real_type>::type_sqrt (t);
    00474 #else
    00475             real_type scale = real_type ();
    00476             real_type sum_squares (1);
    00477             while (-- size >= 0) {
    00478                 real_type u (type_traits<value_type>::norm_2 (*it));
    00479                 if (scale < u) {
    00480                     real_type v (scale / u);
    00481                     sum_squares = sum_squares * v * v + real_type (1);
    00482                     scale = u;
    00483                 } else {
    00484                     real_type v (u / scale);
    00485                     sum_squares += v * v;
    00486                 }
    00487                 ++ it;
    00488             }
    00489             return scale * type_traits<real_type>::type_sqrt (sum_squares);
    00490 #endif
    00491         }
    00492         // Sparse case
    00493         template<class I>
    00494         static BOOST_UBLAS_INLINE
    00495         result_type apply (I it, const I &it_end) {
    00496 #ifndef BOOST_UBLAS_SCALED_NORM
    00497             real_type t = real_type ();
    00498             while (it != it_end) {
    00499                 real_type u (type_traits<value_type>::norm_2 (*it));
    00500                 t +=  u * u;
    00501                 ++ it;
    00502             }
    00503             return type_traits<real_type>::type_sqrt (t);
    00504 #else
    00505             real_type scale = real_type ();
    00506             real_type sum_squares (1);
    00507             while (it != it_end) {
    00508                 real_type u (type_traits<value_type>::norm_2 (*it));
    00509                 if (scale < u) {
    00510                     real_type v (scale / u);
    00511                     sum_squares = sum_squares * v * v + real_type (1);
    00512                     scale = u;
    00513                 } else {
    00514                     real_type v (u / scale);
    00515                     sum_squares += v * v;
    00516                 }
    00517                 ++ it;
    00518             }
    00519             return scale * type_traits<real_type>::type_sqrt (sum_squares);
    00520 #endif
    00521         }
    00522     };
    00523     template<class V>
    00524     struct vector_norm_inf:
    00525         public vector_scalar_real_unary_functor<V> {
    00526         typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
    00527         typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
    00528         typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
    00529
    00530         template<class E>
    00531         static BOOST_UBLAS_INLINE
    00532         result_type apply (const vector_expression<E> &e) {
    00533             real_type t = real_type ();
    00534             typedef typename E::size_type vector_size_type;
    00535             vector_size_type size (e ().size ());
    00536             for (vector_size_type i = 0; i < size; ++ i) {
    00537                 real_type u (type_traits<value_type>::norm_inf (e () (i)));
    00538                 if (u > t)
    00539                     t = u;
    00540             }
    00541             return t;
    00542         }
    00543         // Dense case
    00544         template<class D, class I>
    00545         static BOOST_UBLAS_INLINE
    00546         result_type apply (D size, I it) {
    00547             real_type t = real_type ();
    00548             while (-- size >= 0) {
    00549                 real_type u (type_traits<value_type>::norm_inf (*it));
    00550                 if (u > t)
    00551                     t = u;
    00552                 ++ it;
    00553             }
    00554             return t;
    00555         }
    00556         // Sparse case
    00557         template<class I>
    00558         static BOOST_UBLAS_INLINE
    00559         result_type apply (I it, const I &it_end) {
    00560             real_type t = real_type ();
    00561             while (it != it_end) {
    00562                 real_type u (type_traits<value_type>::norm_inf (*it));
    00563                 if (u > t)
    00564                     t = u;
    00565                 ++ it;
    00566             }
    00567             return t;
    00568         }
    00569     };
    00570
    00571     // Unary returning index
    00572     template<class V>
    00573     struct vector_scalar_index_unary_functor {
    00574         typedef typename V::value_type value_type;
    00575         typedef typename type_traits<value_type>::real_type real_type;
    00576         typedef typename V::size_type result_type;
    00577     };
    00578
    00579     template<class V>
    00580     struct vector_index_norm_inf:
    00581         public vector_scalar_index_unary_functor<V> {
    00582         typedef typename vector_scalar_index_unary_functor<V>::value_type value_type;
    00583         typedef typename vector_scalar_index_unary_functor<V>::real_type real_type;
    00584         typedef typename vector_scalar_index_unary_functor<V>::result_type result_type;
    00585
    00586         template<class E>
    00587         static BOOST_UBLAS_INLINE
    00588         result_type apply (const vector_expression<E> &e) {
    00589             // ISSUE For CBLAS compatibility return 0 index in empty case
    00590             result_type i_norm_inf (0);
    00591             real_type t = real_type ();
    00592             typedef typename E::size_type vector_size_type;
    00593             vector_size_type size (e ().size ());
    00594             for (vector_size_type i = 0; i < size; ++ i) {
    00595                 real_type u (type_traits<value_type>::norm_inf (e () (i)));
    00596                 if (u > t) {
    00597                     i_norm_inf = i;
    00598                     t = u;
    00599                 }
    00600             }
    00601             return i_norm_inf;
    00602         }
    00603         // Dense case
    00604         template<class D, class I>
    00605         static BOOST_UBLAS_INLINE
    00606         result_type apply (D size, I it) {
    00607             // ISSUE For CBLAS compatibility return 0 index in empty case
    00608             result_type i_norm_inf (0);
    00609             real_type t = real_type ();
    00610             while (-- size >= 0) {
    00611                 real_type u (type_traits<value_type>::norm_inf (*it));
    00612                 if (u > t) {
    00613                     i_norm_inf = it.index ();
    00614                     t = u;
    00615                 }
    00616                 ++ it;
    00617             }
    00618             return i_norm_inf;
    00619         }
    00620         // Sparse case
    00621         template<class I>
    00622         static BOOST_UBLAS_INLINE
    00623         result_type apply (I it, const I &it_end) {
    00624             // ISSUE For CBLAS compatibility return 0 index in empty case
    00625             result_type i_norm_inf (0);
    00626             real_type t = real_type ();
    00627             while (it != it_end) {
    00628                 real_type u (type_traits<value_type>::norm_inf (*it));
    00629                 if (u > t) {
    00630                     i_norm_inf = it.index ();
    00631                     t = u;
    00632                 }
    00633                 ++ it;
    00634             }
    00635             return i_norm_inf;
    00636         }
    00637     };
    00638
    00639     // Binary returning scalar
    00640     template<class V1, class V2, class TV>
    00641     struct vector_scalar_binary_functor {
    00642         typedef TV value_type;
    00643         typedef TV result_type;
    00644     };
    00645
    00646     template<class V1, class V2, class TV>
    00647     struct vector_inner_prod:
    00648         public vector_scalar_binary_functor<V1, V2, TV> {
    00649         typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type;
    00650         typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type;
    00651
    00652         template<class C1, class C2>
    00653         static BOOST_UBLAS_INLINE
    00654         result_type apply (const vector_container<C1> &c1,
    00655                            const vector_container<C2> &c2) {
    00656 #ifdef BOOST_UBLAS_USE_SIMD
    00657             using namespace raw;
    00658             typedef typename C1::size_type vector_size_type;
    00659             vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
    00660             const typename V1::value_type *data1 = data_const (c1 ());
    00661             const typename V1::value_type *data2 = data_const (c2 ());
    00662             vector_size_type s1 = stride (c1 ());
    00663             vector_size_type s2 = stride (c2 ());
    00664             result_type t = result_type (0);
    00665             if (s1 == 1 && s2 == 1) {
    00666                 for (vector_size_type i = 0; i < size; ++ i)
    00667                     t += data1 [i] * data2 [i];
    00668             } else if (s2 == 1) {
    00669                 for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
    00670                     t += data1 [i1] * data2 [i];
    00671             } else if (s1 == 1) {
    00672                 for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
    00673                     t += data1 [i] * data2 [i2];
    00674             } else {
    00675                 for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
    00676                     t += data1 [i1] * data2 [i2];
    00677             }
    00678             return t;
    00679 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)
    00680             return boost::numeric::bindings::atlas::dot (c1 (), c2 ());
    00681 #else
    00682             return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2));
    00683 #endif
    00684         }
    00685         template<class E1, class E2>
    00686         static BOOST_UBLAS_INLINE
    00687         result_type apply (const vector_expression<E1> &e1,
    00688                            const vector_expression<E2> &e2) {
    00689             typedef typename E1::size_type vector_size_type;
    00690             vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
    00691             result_type t = result_type (0);
    00692 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
    00693             for (vector_size_type i = 0; i < size; ++ i)
    00694                 t += e1 () (i) * e2 () (i);
    00695 #else
    00696             vector_size_type i (0);
    00697             DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));
    00698 #endif
    00699             return t;
    00700         }
    00701         // Dense case
    00702         template<class D, class I1, class I2>
    00703         static BOOST_UBLAS_INLINE
    00704         result_type apply (D size, I1 it1, I2 it2) {
    00705             result_type t = result_type (0);
    00706 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
    00707             while (-- size >= 0)
    00708                 t += *it1 * *it2, ++ it1, ++ it2;
    00709 #else
    00710             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
    00711 #endif
    00712             return t;
    00713         }
    00714         // Packed case
    00715         template<class I1, class I2>
    00716         static BOOST_UBLAS_INLINE
    00717         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
    00718             result_type t = result_type (0);
    00719             typedef typename I1::difference_type vector_difference_type;
    00720             vector_difference_type it1_size (it1_end - it1);
    00721             vector_difference_type it2_size (it2_end - it2);
    00722             vector_difference_type diff (0);
    00723             if (it1_size > 0 && it2_size > 0)
    00724                 diff = it2.index () - it1.index ();
    00725             if (diff != 0) {
    00726                 vector_difference_type size = (std::min) (diff, it1_size);
    00727                 if (size > 0) {
    00728                     it1 += size;
    00729                     it1_size -= size;
    00730                     diff -= size;
    00731                 }
    00732                 size = (std::min) (- diff, it2_size);
    00733                 if (size > 0) {
    00734                     it2 += size;
    00735                     it2_size -= size;
    00736                     diff += size;
    00737                 }
    00738             }
    00739             vector_difference_type size ((std::min) (it1_size, it2_size));
    00740             while (-- size >= 0)
    00741                 t += *it1 * *it2, ++ it1, ++ it2;
    00742             return t;
    00743         }
    00744         // Sparse case
    00745         template<class I1, class I2>
    00746         static BOOST_UBLAS_INLINE
    00747         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
    00748             result_type t = result_type (0);
    00749             if (it1 != it1_end && it2 != it2_end) {
    00750                 while (true) {
    00751                     if (it1.index () == it2.index ()) {
    00752                         t += *it1 * *it2, ++ it1, ++ it2;
    00753                         if (it1 == it1_end || it2 == it2_end)
    00754                             break;
    00755                     } else if (it1.index () < it2.index ()) {
    00756                         increment (it1, it1_end, it2.index () - it1.index ());
    00757                         if (it1 == it1_end)
    00758                             break;
    00759                     } else if (it1.index () > it2.index ()) {
    00760                         increment (it2, it2_end, it1.index () - it2.index ());
    00761                         if (it2 == it2_end)
    00762                             break;
    00763                     }
    00764                 }
    00765             }
    00766             return t;
    00767         }
    00768     };
    00769
    00770     // Matrix functors
    00771
    00772     // Binary returning vector
    00773     template<class M1, class M2, class TV>
    00774     struct matrix_vector_binary_functor {
    00775         typedef typename M1::size_type size_type;
    00776         typedef typename M1::difference_type difference_type;
    00777         typedef TV value_type;
    00778         typedef TV result_type;
    00779     };
    00780
    00781     template<class M1, class M2, class TV>
    00782     struct matrix_vector_prod1:
    00783         public matrix_vector_binary_functor<M1, M2, TV> {
    00784         typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
    00785         typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
    00786         typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
    00787         typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;
    00788
    00789         template<class C1, class C2>
    00790         static BOOST_UBLAS_INLINE
    00791         result_type apply (const matrix_container<C1> &c1,
    00792                            const vector_container<C2> &c2,
    00793                            size_type i) {
    00794 #ifdef BOOST_UBLAS_USE_SIMD
    00795             using namespace raw;
    00796             size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ());
    00797             const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
    00798             const typename M2::value_type *data2 = data_const (c2 ());
    00799             size_type s1 = stride2 (c1 ());
    00800             size_type s2 = stride (c2 ());
    00801             result_type t = result_type (0);
    00802             if (s1 == 1 && s2 == 1) {
    00803                 for (size_type j = 0; j < size; ++ j)
    00804                     t += data1 [j] * data2 [j];
    00805             } else if (s2 == 1) {
    00806                 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
    00807                     t += data1 [j1] * data2 [j];
    00808             } else if (s1 == 1) {
    00809                 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
    00810                     t += data1 [j] * data2 [j2];
    00811             } else {
    00812                 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
    00813                     t += data1 [j1] * data2 [j2];
    00814             }
    00815             return t;
    00816 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)
    00817             return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ());
    00818 #else
    00819             return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i));
    00820 #endif
    00821         }
    00822         template<class E1, class E2>
    00823         static BOOST_UBLAS_INLINE
    00824         result_type apply (const matrix_expression<E1> &e1,
    00825                            const vector_expression<E2> &e2,
    00826                            size_type i) {
    00827             size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
    00828             result_type t = result_type (0);
    00829 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
    00830             for (size_type j = 0; j < size; ++ j)
    00831                 t += e1 () (i, j) * e2 () (j);
    00832 #else
    00833             size_type j (0);
    00834             DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j));
    00835 #endif
    00836             return t;
    00837         }
    00838         // Dense case
    00839         template<class I1, class I2>
    00840         static BOOST_UBLAS_INLINE
    00841         result_type apply (difference_type size, I1 it1, I2 it2) {
    00842             result_type t = result_type (0);
    00843 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
    00844             while (-- size >= 0)
    00845                 t += *it1 * *it2, ++ it1, ++ it2;
    00846 #else
    00847             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
    00848 #endif
    00849             return t;
    00850         }
    00851         // Packed case
    00852         template<class I1, class I2>
    00853         static BOOST_UBLAS_INLINE
    00854         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
    00855             result_type t = result_type (0);
    00856             difference_type it1_size (it1_end - it1);
    00857             difference_type it2_size (it2_end - it2);
    00858             difference_type diff (0);
    00859             if (it1_size > 0 && it2_size > 0)
    00860                 diff = it2.index () - it1.index2 ();
    00861             if (diff != 0) {
    00862                 difference_type size = (std::min) (diff, it1_size);
    00863                 if (size > 0) {
    00864                     it1 += size;
    00865                     it1_size -= size;
    00866                     diff -= size;
    00867                 }
    00868                 size = (std::min) (- diff, it2_size);
    00869                 if (size > 0) {
    00870                     it2 += size;
    00871                     it2_size -= size;
    00872                     diff += size;
    00873                 }
    00874             }
    00875             difference_type size ((std::min) (it1_size, it2_size));
    00876             while (-- size >= 0)
    00877                 t += *it1 * *it2, ++ it1, ++ it2;
    00878             return t;
    00879         }
    00880         // Sparse case
    00881         template<class I1, class I2>
    00882         static BOOST_UBLAS_INLINE
    00883         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
    00884                            sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
    00885             result_type t = result_type (0);
    00886             if (it1 != it1_end && it2 != it2_end) {
    00887                 size_type it1_index = it1.index2 (), it2_index = it2.index ();
    00888                 while (true) {
    00889                     difference_type compare = it1_index - it2_index;
    00890                     if (compare == 0) {
    00891                         t += *it1 * *it2, ++ it1, ++ it2;
    00892                         if (it1 != it1_end && it2 != it2_end) {
    00893                             it1_index = it1.index2 ();
    00894                             it2_index = it2.index ();
    00895                         } else
    00896                             break;
    00897                     } else if (compare < 0) {
    00898                         increment (it1, it1_end, - compare);
    00899                         if (it1 != it1_end)
    00900                             it1_index = it1.index2 ();
    00901                         else
    00902                             break;
    00903                     } else if (compare > 0) {
    00904                         increment (it2, it2_end, compare);
    00905                         if (it2 != it2_end)
    00906                             it2_index = it2.index ();
    00907                         else
    00908                             break;
    00909                     }
    00910                 }
    00911             }
    00912             return t;
    00913         }
    00914         // Sparse packed case
    00915         template<class I1, class I2>
    00916         static BOOST_UBLAS_INLINE
    00917         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
    00918                            sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
    00919             result_type t = result_type (0);
    00920             while (it1 != it1_end) {
    00921                 t += *it1 * it2 () (it1.index2 ());
    00922                 ++ it1;
    00923             }
    00924             return t;
    00925         }
    00926         // Packed sparse case
    00927         template<class I1, class I2>
    00928         static BOOST_UBLAS_INLINE
    00929         result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
    00930                            packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
    00931             result_type t = result_type (0);
    00932             while (it2 != it2_end) {
    00933                 t += it1 () (it1.index1 (), it2.index ()) * *it2;
    00934                 ++ it2;
    00935             }
    00936             return t;
    00937         }
    00938         // Another dispatcher
    00939         template<class I1, class I2>
    00940         static BOOST_UBLAS_INLINE
    00941         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
    00942                            sparse_bidirectional_iterator_tag) {
    00943             typedef typename I1::iterator_category iterator1_category;
    00944             typedef typename I2::iterator_category iterator2_category;
    00945             return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
    00946         }
    00947     };
    00948
    00949     template<class M1, class M2, class TV>
    00950     struct matrix_vector_prod2:
    00951         public matrix_vector_binary_functor<M1, M2, TV> {
    00952         typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
    00953         typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
    00954         typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
    00955         typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;
    00956
    00957         template<class C1, class C2>
    00958         static BOOST_UBLAS_INLINE
    00959         result_type apply (const vector_container<C1> &c1,
    00960                            const matrix_container<C2> &c2,
    00961                            size_type i) {
    00962 #ifdef BOOST_UBLAS_USE_SIMD
    00963             using namespace raw;
    00964             size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ());
    00965             const typename M1::value_type *data1 = data_const (c1 ());
    00966             const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ());
    00967             size_type s1 = stride (c1 ());
    00968             size_type s2 = stride1 (c2 ());
    00969             result_type t = result_type (0);
    00970             if (s1 == 1 && s2 == 1) {
    00971                 for (size_type j = 0; j < size; ++ j)
    00972                     t += data1 [j] * data2 [j];
    00973             } else if (s2 == 1) {
    00974                 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
    00975                     t += data1 [j1] * data2 [j];
    00976             } else if (s1 == 1) {
    00977                 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
    00978                     t += data1 [j] * data2 [j2];
    00979             } else {
    00980                 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
    00981                     t += data1 [j1] * data2 [j2];
    00982             }
    00983             return t;
    00984 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)
    00985             return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i));
    00986 #else
    00987             return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
    00988 #endif
    00989         }
    00990         template<class E1, class E2>
    00991         static BOOST_UBLAS_INLINE
    00992         result_type apply (const vector_expression<E1> &e1,
    00993                            const matrix_expression<E2> &e2,
    00994                            size_type i) {
    00995             size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
    00996             result_type t = result_type (0);
    00997 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
    00998             for (size_type j = 0; j < size; ++ j)
    00999                 t += e1 () (j) * e2 () (j, i);
    01000 #else
    01001             size_type j (0);
    01002             DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j));
    01003 #endif
    01004             return t;
    01005         }
    01006         // Dense case
    01007         template<class I1, class I2>
    01008         static BOOST_UBLAS_INLINE
    01009         result_type apply (difference_type size, I1 it1, I2 it2) {
    01010             result_type t = result_type (0);
    01011 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
    01012             while (-- size >= 0)
    01013                 t += *it1 * *it2, ++ it1, ++ it2;
    01014 #else
    01015             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
    01016 #endif
    01017             return t;
    01018         }
    01019         // Packed case
    01020         template<class I1, class I2>
    01021         static BOOST_UBLAS_INLINE
    01022         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
    01023             result_type t = result_type (0);
    01024             difference_type it1_size (it1_end - it1);
    01025             difference_type it2_size (it2_end - it2);
    01026             difference_type diff (0);
    01027             if (it1_size > 0 && it2_size > 0)
    01028                 diff = it2.index1 () - it1.index ();
    01029             if (diff != 0) {
    01030                 difference_type size = (std::min) (diff, it1_size);
    01031                 if (size > 0) {
    01032                     it1 += size;
    01033                     it1_size -= size;
    01034                     diff -= size;
    01035                 }
    01036                 size = (std::min) (- diff, it2_size);
    01037                 if (size > 0) {
    01038                     it2 += size;
    01039                     it2_size -= size;
    01040                     diff += size;
    01041                 }
    01042             }
    01043             difference_type size ((std::min) (it1_size, it2_size));
    01044             while (-- size >= 0)
    01045                 t += *it1 * *it2, ++ it1, ++ it2;
    01046             return t;
    01047         }
    01048         // Sparse case
    01049         template<class I1, class I2>
    01050         static BOOST_UBLAS_INLINE
    01051         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
    01052                            sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
    01053             result_type t = result_type (0);
    01054             if (it1 != it1_end && it2 != it2_end) {
    01055                 size_type it1_index = it1.index (), it2_index = it2.index1 ();
    01056                 while (true) {
    01057                     difference_type compare = it1_index - it2_index;
    01058                     if (compare == 0) {
    01059                         t += *it1 * *it2, ++ it1, ++ it2;
    01060                         if (it1 != it1_end && it2 != it2_end) {
    01061                             it1_index = it1.index ();
    01062                             it2_index = it2.index1 ();
    01063                         } else
    01064                             break;
    01065                     } else if (compare < 0) {
    01066                         increment (it1, it1_end, - compare);
    01067                         if (it1 != it1_end)
    01068                             it1_index = it1.index ();
    01069                         else
    01070                             break;
    01071                     } else if (compare > 0) {
    01072                         increment (it2, it2_end, compare);
    01073                         if (it2 != it2_end)
    01074                             it2_index = it2.index1 ();
    01075                         else
    01076                             break;
    01077                     }
    01078                 }
    01079             }
    01080             return t;
    01081         }
    01082         // Packed sparse case
    01083         template<class I1, class I2>
    01084         static BOOST_UBLAS_INLINE
    01085         result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
    01086                            packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
    01087             result_type t = result_type (0);
    01088             while (it2 != it2_end) {
    01089                 t += it1 () (it2.index1 ()) * *it2;
    01090                 ++ it2;
    01091             }
    01092             return t;
    01093         }
    01094         // Sparse packed case
    01095         template<class I1, class I2>
    01096         static BOOST_UBLAS_INLINE
    01097         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
    01098                            sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
    01099             result_type t = result_type (0);
    01100             while (it1 != it1_end) {
    01101                 t += *it1 * it2 () (it1.index (), it2.index2 ());
    01102                 ++ it1;
    01103             }
    01104             return t;
    01105         }
    01106         // Another dispatcher
    01107         template<class I1, class I2>
    01108         static BOOST_UBLAS_INLINE
    01109         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
    01110                            sparse_bidirectional_iterator_tag) {
    01111             typedef typename I1::iterator_category iterator1_category;
    01112             typedef typename I2::iterator_category iterator2_category;
    01113             return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
    01114         }
    01115     };
    01116
    01117     // Binary returning matrix
    01118     template<class M1, class M2, class TV>
    01119     struct matrix_matrix_binary_functor {
    01120         typedef typename M1::size_type size_type;
    01121         typedef typename M1::difference_type difference_type;
    01122         typedef TV value_type;
    01123         typedef TV result_type;
    01124     };
    01125
    01126     template<class M1, class M2, class TV>
    01127     struct matrix_matrix_prod:
    01128         public matrix_matrix_binary_functor<M1, M2, TV> {
    01129         typedef typename matrix_matrix_binary_functor<M1, M2, TV>::size_type size_type;
    01130         typedef typename matrix_matrix_binary_functor<M1, M2, TV>::difference_type difference_type;
    01131         typedef typename matrix_matrix_binary_functor<M1, M2, TV>::value_type value_type;
    01132         typedef typename matrix_matrix_binary_functor<M1, M2, TV>::result_type result_type;
    01133
    01134         template<class C1, class C2>
    01135         static BOOST_UBLAS_INLINE
    01136         result_type apply (const matrix_container<C1> &c1,
    01137                            const matrix_container<C2> &c2,
    01138                            size_type i, size_type j) {
    01139 #ifdef BOOST_UBLAS_USE_SIMD
    01140             using namespace raw;
    01141             size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ());
    01142             const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
    01143             const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ());
    01144             size_type s1 = stride2 (c1 ());
    01145             size_type s2 = stride1 (c2 ());
    01146             result_type t = result_type (0);
    01147             if (s1 == 1 && s2 == 1) {
    01148                 for (size_type k = 0; k < size; ++ k)
    01149                     t += data1 [k] * data2 [k];
    01150             } else if (s2 == 1) {
    01151                 for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1)
    01152                     t += data1 [k1] * data2 [k];
    01153             } else if (s1 == 1) {
    01154                 for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2)
    01155                     t += data1 [k] * data2 [k2];
    01156             } else {
    01157                 for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2)
    01158                     t += data1 [k1] * data2 [k2];
    01159             }
    01160             return t;
    01161 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)
    01162             return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j));
    01163 #else
    01164             return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
    01165 #endif
    01166         }
    01167         template<class E1, class E2>
    01168         static BOOST_UBLAS_INLINE
    01169         result_type apply (const matrix_expression<E1> &e1,
    01170                            const matrix_expression<E2> &e2,
    01171                            size_type i, size_type j) {
    01172             size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
    01173             result_type t = result_type (0);
    01174 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
    01175             for (size_type k = 0; k < size; ++ k)
    01176                 t += e1 () (i, k) * e2 () (k, j);
    01177 #else
    01178             size_type k (0);
    01179             DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k));
    01180 #endif
    01181             return t;
    01182         }
    01183         // Dense case
    01184         template<class I1, class I2>
    01185         static BOOST_UBLAS_INLINE
    01186         result_type apply (difference_type size, I1 it1, I2 it2) {
    01187             result_type t = result_type (0);
    01188 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
    01189             while (-- size >= 0)
    01190                 t += *it1 * *it2, ++ it1, ++ it2;
    01191 #else
    01192             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
    01193 #endif
    01194             return t;
    01195         }
    01196         // Packed case
    01197         template<class I1, class I2>
    01198         static BOOST_UBLAS_INLINE
    01199         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) {
    01200             result_type t = result_type (0);
    01201             difference_type it1_size (it1_end - it1);
    01202             difference_type it2_size (it2_end - it2);
    01203             difference_type diff (0);
    01204             if (it1_size > 0 && it2_size > 0)
    01205                 diff = it2.index1 () - it1.index2 ();
    01206             if (diff != 0) {
    01207                 difference_type size = (std::min) (diff, it1_size);
    01208                 if (size > 0) {
    01209                     it1 += size;
    01210                     it1_size -= size;
    01211                     diff -= size;
    01212                 }
    01213                 size = (std::min) (- diff, it2_size);
    01214                 if (size > 0) {
    01215                     it2 += size;
    01216                     it2_size -= size;
    01217                     diff += size;
    01218                 }
    01219             }
    01220             difference_type size ((std::min) (it1_size, it2_size));
    01221             while (-- size >= 0)
    01222                 t += *it1 * *it2, ++ it1, ++ it2;
    01223             return t;
    01224         }
    01225         // Sparse case
    01226         template<class I1, class I2>
    01227         static BOOST_UBLAS_INLINE
    01228         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
    01229             result_type t = result_type (0);
    01230             if (it1 != it1_end && it2 != it2_end) {
    01231                 size_type it1_index = it1.index2 (), it2_index = it2.index1 ();
    01232                 while (true) {
    01233                     difference_type compare = difference_type (it1_index - it2_index);
    01234                     if (compare == 0) {
    01235                         t += *it1 * *it2, ++ it1, ++ it2;
    01236                         if (it1 != it1_end && it2 != it2_end) {
    01237                             it1_index = it1.index2 ();
    01238                             it2_index = it2.index1 ();
    01239                         } else
    01240                             break;
    01241                     } else if (compare < 0) {
    01242                         increment (it1, it1_end, - compare);
    01243                         if (it1 != it1_end)
    01244                             it1_index = it1.index2 ();
    01245                         else
    01246                             break;
    01247                     } else if (compare > 0) {
    01248                         increment (it2, it2_end, compare);
    01249                         if (it2 != it2_end)
    01250                             it2_index = it2.index1 ();
    01251                         else
    01252                             break;
    01253                     }
    01254                 }
    01255             }
    01256             return t;
    01257         }
    01258     };
    01259
    01260     // Unary returning scalar norm
    01261     template<class M>
    01262     struct matrix_scalar_real_unary_functor {
    01263         typedef typename M::value_type value_type;
    01264         typedef typename type_traits<value_type>::real_type real_type;
    01265         typedef real_type result_type;
    01266     };
    01267
    01268     template<class M>
    01269     struct matrix_norm_1:
    01270         public matrix_scalar_real_unary_functor<M> {
    01271         typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
    01272         typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
    01273         typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
    01274
    01275         template<class E>
    01276         static BOOST_UBLAS_INLINE
    01277         result_type apply (const matrix_expression<E> &e) {
    01278             real_type t = real_type ();
    01279             typedef typename E::size_type matrix_size_type;
    01280             matrix_size_type size2 (e ().size2 ());
    01281             for (matrix_size_type j = 0; j < size2; ++ j) {
    01282                 real_type u = real_type ();
    01283                 matrix_size_type size1 (e ().size1 ());
    01284                 for (matrix_size_type i = 0; i < size1; ++ i) {
    01285                     real_type v (type_traits<value_type>::norm_1 (e () (i, j)));
    01286                     u += v;
    01287                 }
    01288                 if (u > t)
    01289                     t = u;
    01290             }
    01291             return t;
    01292         }
    01293     };
    01294
    01295     template<class M>
    01296     struct matrix_norm_frobenius:
    01297         public matrix_scalar_real_unary_functor<M> {
    01298         typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
    01299         typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
    01300         typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
    01301
    01302         template<class E>
    01303         static BOOST_UBLAS_INLINE
    01304         result_type apply (const matrix_expression<E> &e) {
    01305             real_type t = real_type ();
    01306             typedef typename E::size_type matrix_size_type;
    01307             matrix_size_type size1 (e ().size1 ());
    01308             for (matrix_size_type i = 0; i < size1; ++ i) {
    01309                 matrix_size_type size2 (e ().size2 ());
    01310                 for (matrix_size_type j = 0; j < size2; ++ j) {
    01311                     real_type u (type_traits<value_type>::norm_2 (e () (i, j)));
    01312                     t +=  u * u;
    01313                 }
    01314             }
    01315             return type_traits<real_type>::type_sqrt (t);
    01316         }
    01317     };
    01318
    01319     template<class M>
    01320     struct matrix_norm_inf:
    01321         public matrix_scalar_real_unary_functor<M> {
    01322         typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
    01323         typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
    01324         typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
    01325
    01326         template<class E>
    01327         static BOOST_UBLAS_INLINE
    01328         result_type apply (const matrix_expression<E> &e) {
    01329             real_type t = real_type ();
    01330             typedef typename E::size_type matrix_size_type;
    01331             matrix_size_type size1 (e ().size1 ());
    01332             for (matrix_size_type i = 0; i < size1; ++ i) {
    01333                 real_type u = real_type ();
    01334                 matrix_size_type size2 (e ().size2 ());
    01335                 for (matrix_size_type j = 0; j < size2; ++ j) {
    01336                     real_type v (type_traits<value_type>::norm_inf (e () (i, j)));
    01337                     u += v;
    01338                 }
    01339                 if (u > t)
    01340                     t = u;
    01341             }
    01342             return t;
    01343         }
    01344     };
    01345
    01346     // forward declaration
    01347     template <class Z, class D> struct basic_column_major;
    01348
    01349     // This functor defines storage layout and it's properties
    01350     // matrix (i,j) -> storage [i * size_i + j]
    01351     template <class Z, class D>
    01352     struct basic_row_major {
    01353         typedef Z size_type;
    01354         typedef D difference_type;
    01355         typedef row_major_tag orientation_category;
    01356         typedef basic_column_major<Z,D> transposed_layout;
    01357
    01358         static
    01359         BOOST_UBLAS_INLINE
    01360         size_type storage_size (size_type size_i, size_type size_j) {
    01361             // Guard against size_type overflow
    01362             BOOST_UBLAS_CHECK (size_j == 0 || size_i <= (std::numeric_limits<size_type>::max) () / size_j, bad_size ());
    01363             return size_i * size_j;
    01364         }
    01365
    01366         // Indexing conversion to storage element
    01367         static
    01368         BOOST_UBLAS_INLINE
    01369         size_type element (size_type i, size_type size_i, size_type j, size_type size_j) {
    01370             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
    01371             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
    01372             detail::ignore_unused_variable_warning(size_i);
    01373             // Guard against size_type overflow
    01374             BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ());
    01375             return i * size_j + j;
    01376         }
    01377         static
    01378         BOOST_UBLAS_INLINE
    01379         size_type address (size_type i, size_type size_i, size_type j, size_type size_j) {
    01380             BOOST_UBLAS_CHECK (i <= size_i, bad_index ());
    01381             BOOST_UBLAS_CHECK (j <= size_j, bad_index ());
    01382             // Guard against size_type overflow - address may be size_j past end of storage
    01383             BOOST_UBLAS_CHECK (size_j == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ());
    01384             detail::ignore_unused_variable_warning(size_i);
    01385             return i * size_j + j;
    01386         }
    01387
    01388         // Storage element to index conversion
    01389         static
    01390         BOOST_UBLAS_INLINE
    01391         difference_type distance_i (difference_type k, size_type /* size_i */, size_type size_j) {
    01392             return size_j != 0 ? k / size_j : 0;
    01393         }
    01394         static
    01395         BOOST_UBLAS_INLINE
    01396         difference_type distance_j (difference_type k, size_type /* size_i */, size_type /* size_j */) {
    01397             return k;
    01398         }
    01399         static
    01400         BOOST_UBLAS_INLINE
    01401         size_type index_i (difference_type k, size_type /* size_i */, size_type size_j) {
    01402             return size_j != 0 ? k / size_j : 0;
    01403         }
    01404         static
    01405         BOOST_UBLAS_INLINE
    01406         size_type index_j (difference_type k, size_type /* size_i */, size_type size_j) {
    01407             return size_j != 0 ? k % size_j : 0;
    01408         }
    01409         static
    01410         BOOST_UBLAS_INLINE
    01411         bool fast_i () {
    01412             return false;
    01413         }
    01414         static
    01415         BOOST_UBLAS_INLINE
    01416         bool fast_j () {
    01417             return true;
    01418         }
    01419
    01420         // Iterating storage elements
    01421         template<class I>
    01422         static
    01423         BOOST_UBLAS_INLINE
    01424         void increment_i (I &it, size_type /* size_i */, size_type size_j) {
    01425             it += size_j;
    01426         }
    01427         template<class I>
    01428         static
    01429         BOOST_UBLAS_INLINE
    01430         void increment_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) {
    01431             it += n * size_j;
    01432         }
    01433         template<class I>
    01434         static
    01435         BOOST_UBLAS_INLINE
    01436         void decrement_i (I &it, size_type /* size_i */, size_type size_j) {
    01437             it -= size_j;
    01438         }
    01439         template<class I>
    01440         static
    01441         BOOST_UBLAS_INLINE
    01442         void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) {
    01443             it -= n * size_j;
    01444         }
    01445         template<class I>
    01446         static
    01447         BOOST_UBLAS_INLINE
    01448         void increment_j (I &it, size_type /* size_i */, size_type /* size_j */) {
    01449             ++ it;
    01450         }
    01451         template<class I>
    01452         static
    01453         BOOST_UBLAS_INLINE
    01454         void increment_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
    01455             it += n;
    01456         }
    01457         template<class I>
    01458         static
    01459         BOOST_UBLAS_INLINE
    01460         void decrement_j (I &it, size_type /* size_i */, size_type /* size_j */) {
    01461             -- it;
    01462         }
    01463         template<class I>
    01464         static
    01465         BOOST_UBLAS_INLINE
    01466         void decrement_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
    01467             it -= n;
    01468         }
    01469
    01470         // Triangular access
    01471         static
    01472         BOOST_UBLAS_INLINE
    01473         size_type triangular_size (size_type size_i, size_type size_j) {
    01474             size_type size = (std::max) (size_i, size_j);
    01475             // Guard against size_type overflow - simplified
    01476             BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
    01477             return ((size + 1) * size) / 2;
    01478         }
    01479         static
    01480         BOOST_UBLAS_INLINE
    01481         size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) {
    01482             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
    01483             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
    01484             BOOST_UBLAS_CHECK (i >= j, bad_index ());
    01485             detail::ignore_unused_variable_warning(size_i);
    01486             detail::ignore_unused_variable_warning(size_j);
    01487             // FIXME size_type overflow
    01488             // sigma_i (i + 1) = (i + 1) * i / 2
    01489             // i = 0 1 2 3, sigma = 0 1 3 6
    01490             return ((i + 1) * i) / 2 + j;
    01491         }
    01492         static
    01493         BOOST_UBLAS_INLINE
    01494         size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) {
    01495             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
    01496             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
    01497             BOOST_UBLAS_CHECK (i <= j, bad_index ());
    01498             // FIXME size_type overflow
    01499             // sigma_i (size - i) = size * i - i * (i - 1) / 2
    01500             // i = 0 1 2 3, sigma = 0 4 7 9
    01501             return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i;
    01502         }
    01503
    01504         // Major and minor indices
    01505         static
    01506         BOOST_UBLAS_INLINE
    01507         size_type index_M (size_type index1, size_type /* index2 */) {
    01508             return index1;
    01509         }
    01510         static
    01511         BOOST_UBLAS_INLINE
    01512         size_type index_m (size_type /* index1 */, size_type index2) {
    01513             return index2;
    01514         }
    01515         static
    01516         BOOST_UBLAS_INLINE
    01517         size_type size_M (size_type size_i, size_type /* size_j */) {
    01518             return size_i;
    01519         }
    01520         static
    01521         BOOST_UBLAS_INLINE
    01522         size_type size_m (size_type /* size_i */, size_type size_j) {
    01523             return size_j;
    01524         }
    01525     };
    01526
    01527     // This functor defines storage layout and it's properties
    01528     // matrix (i,j) -> storage [i + j * size_i]
    01529     template <class Z, class D>
    01530     struct basic_column_major {
    01531         typedef Z size_type;
    01532         typedef D difference_type;
    01533         typedef column_major_tag orientation_category;
    01534         typedef basic_row_major<Z,D> transposed_layout;
    01535
    01536         static
    01537         BOOST_UBLAS_INLINE
    01538         size_type storage_size (size_type size_i, size_type size_j) {
    01539             // Guard against size_type overflow
    01540             BOOST_UBLAS_CHECK (size_i == 0 || size_j <= (std::numeric_limits<size_type>::max) () / size_i, bad_size ());
    01541             return size_i * size_j;
    01542         }
    01543
    01544         // Indexing conversion to storage element
    01545         static
    01546         BOOST_UBLAS_INLINE
    01547         size_type element (size_type i, size_type size_i, size_type j, size_type size_j) {
    01548             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
    01549             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
    01550             detail::ignore_unused_variable_warning(size_j);
    01551             // Guard against size_type overflow
    01552             BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ());
    01553             return i + j * size_i;
    01554         }
    01555         static
    01556         BOOST_UBLAS_INLINE
    01557         size_type address (size_type i, size_type size_i, size_type j, size_type size_j) {
    01558             BOOST_UBLAS_CHECK (i <= size_i, bad_index ());
    01559             BOOST_UBLAS_CHECK (j <= size_j, bad_index ());
    01560             detail::ignore_unused_variable_warning(size_j);
    01561             // Guard against size_type overflow - address may be size_i past end of storage
    01562             BOOST_UBLAS_CHECK (size_i == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ());
    01563             return i + j * size_i;
    01564         }
    01565
    01566         // Storage element to index conversion
    01567         static
    01568         BOOST_UBLAS_INLINE
    01569         difference_type distance_i (difference_type k, size_type /* size_i */, size_type /* size_j */) {
    01570             return k;
    01571         }
    01572         static
    01573         BOOST_UBLAS_INLINE
    01574         difference_type distance_j (difference_type k, size_type size_i, size_type /* size_j */) {
    01575             return size_i != 0 ? k / size_i : 0;
    01576         }
    01577         static
    01578         BOOST_UBLAS_INLINE
    01579         size_type index_i (difference_type k, size_type size_i, size_type /* size_j */) {
    01580             return size_i != 0 ? k % size_i : 0;
    01581         }
    01582         static
    01583         BOOST_UBLAS_INLINE
    01584         size_type index_j (difference_type k, size_type size_i, size_type /* size_j */) {
    01585             return size_i != 0 ? k / size_i : 0;
    01586         }
    01587         static
    01588         BOOST_UBLAS_INLINE
    01589         bool fast_i () {
    01590             return true;
    01591         }
    01592         static
    01593         BOOST_UBLAS_INLINE
    01594         bool fast_j () {
    01595             return false;
    01596         }
    01597
    01598         // Iterating
    01599         template<class I>
    01600         static
    01601         BOOST_UBLAS_INLINE
    01602         void increment_i (I &it, size_type /* size_i */, size_type /* size_j */) {
    01603             ++ it;
    01604         }
    01605         template<class I>
    01606         static
    01607         BOOST_UBLAS_INLINE
    01608         void increment_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
    01609             it += n;
    01610         }
    01611         template<class I>
    01612         static
    01613         BOOST_UBLAS_INLINE
    01614         void decrement_i (I &it, size_type /* size_i */, size_type /* size_j */) {
    01615             -- it;
    01616         }
    01617         template<class I>
    01618         static
    01619         BOOST_UBLAS_INLINE
    01620         void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
    01621             it -= n;
    01622         }
    01623         template<class I>
    01624         static
    01625         BOOST_UBLAS_INLINE
    01626         void increment_j (I &it, size_type size_i, size_type /* size_j */) {
    01627             it += size_i;
    01628         }
    01629         template<class I>
    01630         static
    01631         BOOST_UBLAS_INLINE
    01632         void increment_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) {
    01633             it += n * size_i;
    01634         }
    01635         template<class I>
    01636         static
    01637         BOOST_UBLAS_INLINE
    01638         void decrement_j (I &it, size_type size_i, size_type /* size_j */) {
    01639             it -= size_i;
    01640         }
    01641         template<class I>
    01642         static
    01643         BOOST_UBLAS_INLINE
    01644         void decrement_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) {
    01645             it -= n* size_i;
    01646         }
    01647
    01648         // Triangular access
    01649         static
    01650         BOOST_UBLAS_INLINE
    01651         size_type triangular_size (size_type size_i, size_type size_j) {
    01652             size_type size = (std::max) (size_i, size_j);
    01653             // Guard against size_type overflow - simplified
    01654             BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
    01655             return ((size + 1) * size) / 2;
    01656         }
    01657         static
    01658         BOOST_UBLAS_INLINE
    01659         size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) {
    01660             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
    01661             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
    01662             BOOST_UBLAS_CHECK (i >= j, bad_index ());
    01663             // FIXME size_type overflow
    01664             // sigma_j (size - j) = size * j - j * (j - 1) / 2
    01665             // j = 0 1 2 3, sigma = 0 4 7 9
    01666             return i - j + (j * (2 * (std::max) (size_i, size_j) - j + 1)) / 2;
    01667         }
    01668         static
    01669         BOOST_UBLAS_INLINE
    01670         size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) {
    01671             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
    01672             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
    01673             BOOST_UBLAS_CHECK (i <= j, bad_index ());
    01674             // FIXME size_type overflow
    01675             // sigma_j (j + 1) = (j + 1) * j / 2
    01676             // j = 0 1 2 3, sigma = 0 1 3 6
    01677             return i + ((j + 1) * j) / 2;
    01678         }
    01679
    01680         // Major and minor indices
    01681         static
    01682         BOOST_UBLAS_INLINE
    01683         size_type index_M (size_type /* index1 */, size_type index2) {
    01684             return index2;
    01685         }
    01686         static
    01687         BOOST_UBLAS_INLINE
    01688         size_type index_m (size_type index1, size_type /* index2 */) {
    01689             return index1;
    01690         }
    01691         static
    01692         BOOST_UBLAS_INLINE
    01693         size_type size_M (size_type /* size_i */, size_type size_j) {
    01694             return size_j;
    01695         }
    01696         static
    01697         BOOST_UBLAS_INLINE
    01698         size_type size_m (size_type size_i, size_type /* size_j */) {
    01699             return size_i;
    01700         }
    01701     };
    01702
    01703
    01704     template <class Z>
    01705     struct basic_full {
    01706         typedef Z size_type;
    01707
    01708         template<class L>
    01709         static
    01710         BOOST_UBLAS_INLINE
    01711         size_type packed_size (L, size_type size_i, size_type size_j) {
    01712             return L::storage_size (size_i, size_j);
    01713         }
    01714
    01715         static
    01716         BOOST_UBLAS_INLINE
    01717         bool zero (size_type /* i */, size_type /* j */) {
    01718             return false;
    01719         }
    01720         static
    01721         BOOST_UBLAS_INLINE
    01722         bool one (size_type /* i */, size_type /* j */) {
    01723             return false;
    01724         }
    01725         static
    01726         BOOST_UBLAS_INLINE
    01727         bool other (size_type /* i */, size_type /* j */) {
    01728             return true;
    01729         }
    01730         // FIXME: this should not be used at all
    01731         static
    01732         BOOST_UBLAS_INLINE
    01733         size_type restrict1 (size_type i, size_type /* j */) {
    01734             return i;
    01735         }
    01736         static
    01737         BOOST_UBLAS_INLINE
    01738         size_type restrict2 (size_type /* i */, size_type j) {
    01739             return j;
    01740         }
    01741         static
    01742         BOOST_UBLAS_INLINE
    01743         size_type mutable_restrict1 (size_type i, size_type /* j */) {
    01744             return i;
    01745         }
    01746         static
    01747         BOOST_UBLAS_INLINE
    01748         size_type mutable_restrict2 (size_type /* i */, size_type j) {
    01749             return j;
    01750         }
    01751     };
    01752
    01753     namespace detail {
    01754         template < class L >
    01755         struct transposed_structure {
    01756             typedef typename L::size_type size_type;
    01757
    01758             template<class LAYOUT>
    01759             static
    01760             BOOST_UBLAS_INLINE
    01761             size_type packed_size (LAYOUT l, size_type size_i, size_type size_j) {
    01762                 return L::packed_size(l, size_j, size_i);
    01763             }
    01764
    01765             static
    01766             BOOST_UBLAS_INLINE
    01767             bool zero (size_type i, size_type j) {
    01768                 return L::zero(j, i);
    01769             }
    01770             static
    01771             BOOST_UBLAS_INLINE
    01772             bool one (size_type i, size_type j) {
    01773                 return L::one(j, i);
    01774             }
    01775             static
    01776             BOOST_UBLAS_INLINE
    01777             bool other (size_type i, size_type j) {
    01778                 return L::other(j, i);
    01779             }
    01780             template<class LAYOUT>
    01781             static
    01782             BOOST_UBLAS_INLINE
    01783             size_type element (LAYOUT /* l */, size_type i, size_type size_i, size_type j, size_type size_j) {
    01784                 return L::element(typename LAYOUT::transposed_layout(), j, size_j, i, size_i);
    01785             }
    01786
    01787             static
    01788             BOOST_UBLAS_INLINE
    01789             size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
    01790                 return L::restrict2(j, i, size2, size1);
    01791             }
    01792             static
    01793             BOOST_UBLAS_INLINE
    01794             size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
    01795                 return L::restrict1(j, i, size2, size1);
    01796             }
    01797             static
    01798             BOOST_UBLAS_INLINE
    01799             size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
    01800                 return L::mutable_restrict2(j, i, size2, size1);
    01801             }
    01802             static
    01803             BOOST_UBLAS_INLINE
    01804             size_type mutable_restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
    01805                 return L::mutable_restrict1(j, i, size2, size1);
    01806             }
    01807
    01808             static
    01809             BOOST_UBLAS_INLINE
    01810             size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
    01811                 return L::global_restrict2(index2, size2, index1, size1);
    01812             }
    01813             static
    01814             BOOST_UBLAS_INLINE
    01815             size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
    01816                 return L::global_restrict1(index2, size2, index1, size1);
    01817             }
    01818             static
    01819             BOOST_UBLAS_INLINE
    01820             size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
    01821                 return L::global_mutable_restrict2(index2, size2, index1, size1);
    01822             }
    01823             static
    01824             BOOST_UBLAS_INLINE
    01825             size_type global_mutable_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
    01826                 return L::global_mutable_restrict1(index2, size2, index1, size1);
    01827             }
    01828         };
    01829     }
    01830
    01831     template <class Z>
    01832     struct basic_lower {
    01833         typedef Z size_type;
    01834         typedef lower_tag triangular_type;
    01835
    01836         template<class L>
    01837         static
    01838         BOOST_UBLAS_INLINE
    01839         size_type packed_size (L, size_type size_i, size_type size_j) {
    01840             return L::triangular_size (size_i, size_j);
    01841         }
    01842
    01843         static
    01844         BOOST_UBLAS_INLINE
    01845         bool zero (size_type i, size_type j) {
    01846             return j > i;
    01847         }
    01848         static
    01849         BOOST_UBLAS_INLINE
    01850         bool one (size_type /* i */, size_type /* j */) {
    01851             return false;
    01852         }
    01853         static
    01854         BOOST_UBLAS_INLINE
    01855         bool other (size_type i, size_type j) {
    01856             return j <= i;
    01857         }
    01858         template<class L>
    01859         static
    01860         BOOST_UBLAS_INLINE
    01861         size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
    01862             return L::lower_element (i, size_i, j, size_j);
    01863         }
    01864
    01865         // return nearest valid index in column j
    01866         static
    01867         BOOST_UBLAS_INLINE
    01868         size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
    01869             return (std::max)(j, (std::min) (size1, i));
    01870         }
    01871         // return nearest valid index in row i
    01872         static
    01873         BOOST_UBLAS_INLINE
    01874         size_type restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
    01875             return (std::max)(size_type(0), (std::min) (i+1, j));
    01876         }
    01877         // return nearest valid mutable index in column j
    01878         static
    01879         BOOST_UBLAS_INLINE
    01880         size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
    01881             return (std::max)(j, (std::min) (size1, i));
    01882         }
    01883         // return nearest valid mutable index in row i
    01884         static
    01885         BOOST_UBLAS_INLINE
    01886         size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
    01887             return (std::max)(size_type(0), (std::min) (i+1, j));
    01888         }
    01889
    01890         // return an index between the first and (1+last) filled row
    01891         static
    01892         BOOST_UBLAS_INLINE
    01893         size_type global_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
    01894             return (std::max)(size_type(0), (std::min)(size1, index1) );
    01895         }
    01896         // return an index between the first and (1+last) filled column
    01897         static
    01898         BOOST_UBLAS_INLINE
    01899         size_type global_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
    01900             return (std::max)(size_type(0), (std::min)(size2, index2) );
    01901         }
    01902
    01903         // return an index between the first and (1+last) filled mutable row
    01904         static
    01905         BOOST_UBLAS_INLINE
    01906         size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
    01907             return (std::max)(size_type(0), (std::min)(size1, index1) );
    01908         }
    01909         // return an index between the first and (1+last) filled mutable column
    01910         static
    01911         BOOST_UBLAS_INLINE
    01912         size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
    01913             return (std::max)(size_type(0), (std::min)(size2, index2) );
    01914         }
    01915     };
    01916
    01917     // the first row only contains a single 1. Thus it is not stored.
    01918     template <class Z>
    01919     struct basic_unit_lower : public basic_lower<Z> {
    01920         typedef Z size_type;
    01921         typedef unit_lower_tag triangular_type;
    01922
    01923         template<class L>
    01924         static
    01925         BOOST_UBLAS_INLINE
    01926         size_type packed_size (L, size_type size_i, size_type size_j) {
    01927             // Zero size strict triangles are bad at this point
    01928             BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ());
    01929             return L::triangular_size (size_i - 1, size_j - 1);
    01930         }
    01931
    01932         static
    01933         BOOST_UBLAS_INLINE
    01934         bool one (size_type i, size_type j) {
    01935             return j == i;
    01936         }
    01937         static
    01938         BOOST_UBLAS_INLINE
    01939         bool other (size_type i, size_type j) {
    01940             return j < i;
    01941         }
    01942         template<class L>
    01943         static
    01944         BOOST_UBLAS_INLINE
    01945         size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
    01946             // Zero size strict triangles are bad at this point
    01947             BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ());
    01948             return L::lower_element (i-1, size_i - 1, j, size_j - 1);
    01949         }
    01950
    01951         static
    01952         BOOST_UBLAS_INLINE
    01953         size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
    01954             return (std::max)(j+1, (std::min) (size1, i));
    01955         }
    01956         static
    01957         BOOST_UBLAS_INLINE
    01958         size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
    01959             return (std::max)(size_type(0), (std::min) (i, j));
    01960         }
    01961
    01962         // return an index between the first and (1+last) filled mutable row
    01963         static
    01964         BOOST_UBLAS_INLINE
    01965         size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
    01966             return (std::max)(size_type(1), (std::min)(size1, index1) );
    01967         }
    01968         // return an index between the first and (1+last) filled mutable column
    01969         static
    01970         BOOST_UBLAS_INLINE
    01971         size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
    01972             BOOST_UBLAS_CHECK( size2 >= 1 , external_logic() );
    01973             return (std::max)(size_type(0), (std::min)(size2-1, index2) );
    01974         }
    01975     };
    01976
    01977     // the first row only contains no element. Thus it is not stored.
    01978     template <class Z>
    01979     struct basic_strict_lower : public basic_unit_lower<Z> {
    01980         typedef Z size_type;
    01981         typedef strict_lower_tag triangular_type;
    01982
    01983         template<class L>
    01984         static
    01985         BOOST_UBLAS_INLINE
    01986         size_type packed_size (L, size_type size_i, size_type size_j) {
    01987             // Zero size strict triangles are bad at this point
    01988             BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ());
    01989             return L::triangular_size (size_i - 1, size_j - 1);
    01990         }
    01991
    01992         static
    01993         BOOST_UBLAS_INLINE
    01994         bool zero (size_type i, size_type j) {
    01995             return j >= i;
    01996         }
    01997         static
    01998         BOOST_UBLAS_INLINE
    01999         bool one (size_type /*i*/, size_type /*j*/) {
    02000             return false;
    02001         }
    02002         static
    02003         BOOST_UBLAS_INLINE
    02004         bool other (size_type i, size_type j) {
    02005             return j < i;
    02006         }
    02007         template<class L>
    02008         static
    02009         BOOST_UBLAS_INLINE
    02010         size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
    02011             // Zero size strict triangles are bad at this point
    02012             BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ());
    02013             return L::lower_element (i-1, size_i - 1, j, size_j - 1);
    02014         }
    02015
    02016         static
    02017         BOOST_UBLAS_INLINE
    02018         size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
    02019             return basic_unit_lower<Z>::mutable_restrict1(i, j, size1, size2);
    02020         }
    02021         static
    02022         BOOST_UBLAS_INLINE
    02023         size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
    02024             return basic_unit_lower<Z>::mutable_restrict2(i, j, size1, size2);
    02025         }
    02026
    02027         // return an index between the first and (1+last) filled row
    02028         static
    02029         BOOST_UBLAS_INLINE
    02030         size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
    02031             return basic_unit_lower<Z>::global_mutable_restrict1(index1, size1, index2, size2);
    02032         }
    02033         // return an index between the first and (1+last) filled column
    02034         static
    02035         BOOST_UBLAS_INLINE
    02036         size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
    02037             return basic_unit_lower<Z>::global_mutable_restrict2(index1, size1, index2, size2);
    02038         }
    02039     };
    02040
    02041
    02042     template <class Z>
    02043     struct basic_upper : public detail::transposed_structure<basic_lower<Z> >
    02044     {
    02045         typedef upper_tag triangular_type;
    02046     };
    02047
    02048     template <class Z>
    02049     struct basic_unit_upper : public detail::transposed_structure<basic_unit_lower<Z> >
    02050     {
    02051         typedef unit_upper_tag triangular_type;
    02052     };
    02053
    02054     template <class Z>
    02055     struct basic_strict_upper : public detail::transposed_structure<basic_strict_lower<Z> >
    02056     {
    02057         typedef strict_upper_tag triangular_type;
    02058     };
    02059
    02060
    02061 }}}
    02062
    02063 #endif