Boost.uBlas 1.49
Linear Algebra in C++: matrices, vectors and numeric algorithms

functional.hpp

Go to the documentation of this file.
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