Symaptic: os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/functional.hpp@260cb5ec6c19

     1 //

     2 //  Copyright (c) 2000-2002

     3 //  Joerg Walter, Mathias Koch

     4 //

     5 //  Permission to use, copy, modify, distribute and sell this software

     6 //  and its documentation for any purpose is hereby granted without fee,

     7 //  provided that the above copyright notice appear in all copies and

     8 //  that both that copyright notice and this permission notice appear

     9 //  in supporting documentation.  The authors make no representations

    10 //  about the suitability of this software for any purpose.

    11 //  It is provided "as is" without express or implied warranty.

    12 //

    13 //  The authors gratefully acknowledge the support of

    14 //  GeNeSys mbH & Co. KG in producing this work.

    15 //

    17 #ifndef _BOOST_UBLAS_FUNCTIONAL_

    18 #define _BOOST_UBLAS_FUNCTIONAL_

    20 #include <functional>

    22 #include <boost/numeric/ublas/traits.hpp>

    23 #ifdef BOOST_UBLAS_USE_DUFF_DEVICE

    24 #include <boost/numeric/ublas/detail/duff.hpp>

    25 #endif

    26 #ifdef BOOST_UBLAS_USE_SIMD

    27 #include <boost/numeric/ublas/detail/raw.hpp>

    28 #else

    29 namespace boost { namespace numeric { namespace ublas { namespace raw {

    30 }}}}

    31 #endif

    32 #ifdef BOOST_UBLAS_HAVE_BINDINGS

    33 #include <boost/numeric/bindings/traits/std_vector.hpp>

    34 #include <boost/numeric/bindings/traits/ublas_vector.hpp>

    35 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>

    36 #include <boost/numeric/bindings/atlas/cblas.hpp>

    37 #endif

    39 #include <boost/numeric/ublas/detail/definitions.hpp>

    43 namespace boost { namespace numeric { namespace ublas {

    45     // Scalar functors

    47     // Unary

    48     template<class T>

    49     struct scalar_unary_functor {

    50         typedef T value_type;

    51         typedef typename type_traits<T>::const_reference argument_type;

    52         typedef typename type_traits<T>::value_type result_type;

    53     };

    55     template<class T>

    56     struct scalar_identity:

    57         public scalar_unary_functor<T> {

    58         typedef typename scalar_unary_functor<T>::argument_type argument_type;

    59         typedef typename scalar_unary_functor<T>::result_type result_type;

    61         static BOOST_UBLAS_INLINE

    62         result_type apply (argument_type t) {

    63             return t;

    64         }

    65     };

    66     template<class T>

    67     struct scalar_negate:

    68         public scalar_unary_functor<T> {

    69         typedef typename scalar_unary_functor<T>::argument_type argument_type;

    70         typedef typename scalar_unary_functor<T>::result_type result_type;

    72         static BOOST_UBLAS_INLINE

    73         result_type apply (argument_type t) {

    74             return - t;

    75         }

    76     };

    77     template<class T>

    78     struct scalar_conj:

    79         public scalar_unary_functor<T> {

    80         typedef typename scalar_unary_functor<T>::value_type value_type;

    81         typedef typename scalar_unary_functor<T>::argument_type argument_type;

    82         typedef typename scalar_unary_functor<T>::result_type result_type;

    84         static BOOST_UBLAS_INLINE

    85         result_type apply (argument_type t) {

    86             return type_traits<value_type>::conj (t);

    87         }

    88     };

    90     // Unary returning real

    91     template<class T>

    92     struct scalar_real_unary_functor {

    93         typedef T value_type;

    94         typedef typename type_traits<T>::const_reference argument_type;

    95         typedef typename type_traits<T>::real_type result_type;

    96     };

    98     template<class T>

    99     struct scalar_real:

   100         public scalar_real_unary_functor<T> {

   101         typedef typename scalar_real_unary_functor<T>::value_type value_type;

   102         typedef typename scalar_real_unary_functor<T>::argument_type argument_type;

   103         typedef typename scalar_real_unary_functor<T>::result_type result_type;

   105         static BOOST_UBLAS_INLINE

   106         result_type apply (argument_type t) {

   107             return type_traits<value_type>::real (t);

   108         }

   109     };

   110     template<class T>

   111     struct scalar_imag:

   112         public scalar_real_unary_functor<T> {

   113         typedef typename scalar_real_unary_functor<T>::value_type value_type;

   114         typedef typename scalar_real_unary_functor<T>::argument_type argument_type;

   115         typedef typename scalar_real_unary_functor<T>::result_type result_type;

   117         static BOOST_UBLAS_INLINE

   118         result_type apply (argument_type t) {

   119             return type_traits<value_type>::imag (t);

   120         }

   121     };

   123     // Binary

   124     template<class T1, class T2>

   125     struct scalar_binary_functor {

   126         typedef typename type_traits<T1>::const_reference argument1_type;

   127         typedef typename type_traits<T2>::const_reference argument2_type;

   128         typedef typename promote_traits<T1, T2>::promote_type result_type;

   129     };

   131     template<class T1, class T2>

   132     struct scalar_plus:

   133         public scalar_binary_functor<T1, T2> {

   134         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;

   135         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;

   136         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;

   138         static BOOST_UBLAS_INLINE

   139         result_type apply (argument1_type t1, argument2_type t2) {

   140             return t1 + t2;

   141         }

   142     };

   143     template<class T1, class T2>

   144     struct scalar_minus:

   145         public scalar_binary_functor<T1, T2> {

   146         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;

   147         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;

   148         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;

   150         static BOOST_UBLAS_INLINE

   151         result_type apply (argument1_type t1, argument2_type t2) {

   152             return t1 - t2;

   153         }

   154     };

   155     template<class T1, class T2>

   156     struct scalar_multiplies:

   157         public scalar_binary_functor<T1, T2> {

   158         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;

   159         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;

   160         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;

   162         static BOOST_UBLAS_INLINE

   163         result_type apply (argument1_type t1, argument2_type t2) {

   164             return t1 * t2;

   165         }

   166     };

   167     template<class T1, class T2>

   168     struct scalar_divides:

   169         public scalar_binary_functor<T1, T2> {

   170         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;

   171         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;

   172         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;

   174         static BOOST_UBLAS_INLINE

   175         result_type apply (argument1_type t1, argument2_type t2) {

   176             return t1 / t2;

   177         }

   178     };

   180     template<class T1, class T2>

   181     struct scalar_binary_assign_functor {

   182         // ISSUE Remove reference to avoid reference to reference problems

   183         typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;

   184         typedef typename type_traits<T2>::const_reference argument2_type;

   185     };

   187     struct assign_tag {};

   188     struct computed_assign_tag {};

   190     template<class T1, class T2>

   191     struct scalar_assign:

   192         public scalar_binary_assign_functor<T1, T2> {

   193         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;

   194         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;

   195 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )

   196         static const bool computed ;

   197 #else

   198         static const bool computed = false ;

   199 #endif

   201         static BOOST_UBLAS_INLINE

   202         void apply (argument1_type t1, argument2_type t2) {

   203             t1 = t2;

   204         }

   206         template<class U1, class U2>

   207         struct rebind {

   208             typedef scalar_assign<U1, U2> other;

   209         };

   210     };

   212 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )

   213     template<class T1, class T2>

   214     const bool scalar_assign<T1,T2>::computed = false;

   215 #endif

   217     template<class T1, class T2>

   218     struct scalar_plus_assign:

   219         public scalar_binary_assign_functor<T1, T2> {

   220         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;

   221         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;

   222 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )

   223         static const bool computed ;

   224 #else

   225       static const bool computed = true ;

   226 #endif

   228         static BOOST_UBLAS_INLINE

   229         void apply (argument1_type t1, argument2_type t2) {

   230             t1 += t2;

   231         }

   233         template<class U1, class U2>

   234         struct rebind {

   235             typedef scalar_plus_assign<U1, U2> other;

   236         };

   237     };

   239 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )

   240     template<class T1, class T2>

   241     const bool scalar_plus_assign<T1,T2>::computed = true;

   242 #endif

   244     template<class T1, class T2>

   245     struct scalar_minus_assign:

   246         public scalar_binary_assign_functor<T1, T2> {

   247         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;

   248         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;

   249 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )

   250         static const bool computed ;

   251 #else

   252         static const bool computed = true ;

   253 #endif

   255         static BOOST_UBLAS_INLINE

   256         void apply (argument1_type t1, argument2_type t2) {

   257             t1 -= t2;

   258         }

   260         template<class U1, class U2>

   261         struct rebind {

   262             typedef scalar_minus_assign<U1, U2> other;

   263         };

   264     };

   266 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )

   267     template<class T1, class T2>

   268     const bool scalar_minus_assign<T1,T2>::computed = true;

   269 #endif

   271     template<class T1, class T2>

   272     struct scalar_multiplies_assign:

   273         public scalar_binary_assign_functor<T1, T2> {

   274         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;

   275         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;

   276         static const bool computed = true;

   278         static BOOST_UBLAS_INLINE

   279         void apply (argument1_type t1, argument2_type t2) {

   280             t1 *= t2;

   281         }

   283         template<class U1, class U2>

   284         struct rebind {

   285             typedef scalar_multiplies_assign<U1, U2> other;

   286         };

   287     };

   288     template<class T1, class T2>

   289     struct scalar_divides_assign:

   290         public scalar_binary_assign_functor<T1, T2> {

   291         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;

   292         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;

   293         static const bool computed ;

   295         static BOOST_UBLAS_INLINE

   296         void apply (argument1_type t1, argument2_type t2) {

   297             t1 /= t2;

   298         }

   300         template<class U1, class U2>

   301         struct rebind {

   302             typedef scalar_divides_assign<U1, U2> other;

   303         };

   304     };

   305     template<class T1, class T2>

   306     const bool scalar_divides_assign<T1,T2>::computed = true;

   308     template<class T1, class T2>

   309     struct scalar_binary_swap_functor {

   310         typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;

   311         typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type;

   312     };

   314     template<class T1, class T2>

   315     struct scalar_swap:

   316         public scalar_binary_swap_functor<T1, T2> {

   317         typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type;

   318         typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type;

   320         static BOOST_UBLAS_INLINE

   321         void apply (argument1_type t1, argument2_type t2) {

   322             std::swap (t1, t2);

   323         }

   325         template<class U1, class U2>

   326         struct rebind {

   327             typedef scalar_swap<U1, U2> other;

   328         };

   329     };

   331     // Vector functors

   333     // Unary returning scalar

   334     template<class T>

   335     struct vector_scalar_unary_functor {

   336         typedef std::size_t size_type;

   337         typedef std::ptrdiff_t difference_type;

   338         typedef T value_type;

   339         typedef T result_type;

   340     };

   342     template<class T>

   343     struct vector_sum:

   344         public vector_scalar_unary_functor<T> {

   345         typedef typename vector_scalar_unary_functor<T>::size_type size_type;

   346         typedef typename vector_scalar_unary_functor<T>::difference_type difference_type;

   347         typedef typename vector_scalar_unary_functor<T>::value_type value_type;

   348         typedef typename vector_scalar_unary_functor<T>::result_type result_type;

   350         template<class E>

   351         static BOOST_UBLAS_INLINE

   352         result_type apply (const vector_expression<E> &e) {

   353             result_type t = result_type (0);

   354             size_type size (e ().size ());

   355             for (size_type i = 0; i < size; ++ i)

   356                 t += e () (i);

   357             return t;

   358         }

   359         // Dense case

   360         template<class I>

   361         static BOOST_UBLAS_INLINE

   362         result_type apply (difference_type size, I it) {

   363             result_type t = result_type (0);

   364             while (-- size >= 0)

   365                 t += *it, ++ it;

   366             return t;

   367         }

   368         // Sparse case

   369         template<class I>

   370         static BOOST_UBLAS_INLINE

   371         result_type apply (I it, const I &it_end) {

   372             result_type t = result_type (0);

   373             while (it != it_end)

   374                 t += *it, ++ it;

   375             return t;

   376         }

   377     };

   379     // Unary returning real scalar

   380     template<class T>

   381     struct vector_scalar_real_unary_functor {

   382         typedef std::size_t size_type;

   383         typedef std::ptrdiff_t difference_type;

   384         typedef T value_type;

   385         typedef typename type_traits<T>::real_type real_type;

   386         typedef real_type result_type;

   387     };

   389     template<class T>

   390     struct vector_norm_1:

   391         public vector_scalar_real_unary_functor<T> {

   392         typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;

   393         typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;

   394         typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;

   395         typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;

   396         typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;

   398         template<class E>

   399         static BOOST_UBLAS_INLINE

   400         result_type apply (const vector_expression<E> &e) {

   401             real_type t = real_type ();

   402             size_type size (e ().size ());

   403             for (size_type i = 0; i < size; ++ i) {

   404                 real_type u (type_traits<value_type>::type_abs (e () (i)));

   405                 t += u;

   406             }

   407             return t;

   408         }

   409         // Dense case

   410         template<class I>

   411         static BOOST_UBLAS_INLINE

   412         result_type apply (difference_type size, I it) {

   413             real_type t = real_type ();

   414             while (-- size >= 0) {

   415                 real_type u (type_traits<value_type>::norm_1 (*it));

   416                 t += u;

   417                 ++ it;

   418             }

   419             return t;

   420         }

   421         // Sparse case

   422         template<class I>

   423         static BOOST_UBLAS_INLINE

   424         result_type apply (I it, const I &it_end) {

   425             real_type t = real_type ();

   426             while (it != it_end) {

   427                 real_type u (type_traits<value_type>::norm_1 (*it));

   428                 t += u;

   429                 ++ it;

   430             }

   431             return t;

   432         }

   433     };

   434     template<class T>

   435     struct vector_norm_2:

   436         public vector_scalar_real_unary_functor<T> {

   437         typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;

   438         typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;

   439         typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;

   440         typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;

   441         typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;

   443         template<class E>

   444         static BOOST_UBLAS_INLINE

   445         result_type apply (const vector_expression<E> &e) {

   446 #ifndef BOOST_UBLAS_SCALED_NORM

   447             real_type t = real_type ();

   448             size_type size (e ().size ());

   449             for (size_type i = 0; i < size; ++ i) {

   450                 real_type u (type_traits<value_type>::norm_2 (e () (i)));

   451                 t +=  u * u;

   452             }

   453             return type_traits<real_type>::type_sqrt (t);

   454 #else

   455             real_type scale = real_type ();

   456             real_type sum_squares (1);

   457             size_type size (e ().size ());

   458             for (size_type i = 0; i < size; ++ i) {

   459                 real_type u (type_traits<value_type>::norm_2 (e () (i)));

   460                 if (scale < u) {

   461                     real_type v (scale / u);

   462                     sum_squares = sum_squares * v * v + real_type (1);

   463                     scale = u;

   464                 } else {

   465                     real_type v (u / scale);

   466                     sum_squares += v * v;

   467                 }

   468             }

   469             return scale * type_traits<real_type>::type_sqrt (sum_squares);

   470 #endif

   471         }

   472         // Dense case

   473         template<class I>

   474         static BOOST_UBLAS_INLINE

   475         result_type apply (difference_type size, I it) {

   476 #ifndef BOOST_UBLAS_SCALED_NORM

   477             real_type t = real_type ();

   478             while (-- size >= 0) {

   479                 real_type u (type_traits<value_type>::norm_2 (*it));

   480                 t +=  u * u;

   481                 ++ it;

   482             }

   483             return type_traits<real_type>::type_sqrt (t);

   484 #else

   485             real_type scale = real_type ();

   486             real_type sum_squares (1);

   487             while (-- size >= 0) {

   488                 real_type u (type_traits<value_type>::norm_2 (*it));

   489                 if (scale < u) {

   490                     real_type v (scale / u);

   491                     sum_squares = sum_squares * v * v + real_type (1);

   492                     scale = u;

   493                 } else {

   494                     real_type v (u / scale);

   495                     sum_squares += v * v;

   496                 }

   497                 ++ it;

   498             }

   499             return scale * type_traits<real_type>::type_sqrt (sum_squares);

   500 #endif

   501         }

   502         // Sparse case

   503         template<class I>

   504         static BOOST_UBLAS_INLINE

   505         result_type apply (I it, const I &it_end) {

   506 #ifndef BOOST_UBLAS_SCALED_NORM

   507             real_type t = real_type ();

   508             while (it != it_end) {

   509                 real_type u (type_traits<value_type>::norm_2 (*it));

   510                 t +=  u * u;

   511                 ++ it;

   512             }

   513             return type_traits<real_type>::type_sqrt (t);

   514 #else

   515             real_type scale = real_type ();

   516             real_type sum_squares (1);

   517             while (it != it_end) {

   518                 real_type u (type_traits<value_type>::norm_2 (*it));

   519                 if (scale < u) {

   520                     real_type v (scale / u);

   521                     sum_squares = sum_squares * v * v + real_type (1);

   522                     scale = u;

   523                 } else {

   524                     real_type v (u / scale);

   525                     sum_squares += v * v;

   526                 }

   527                 ++ it;

   528             }

   529             return scale * type_traits<real_type>::type_sqrt (sum_squares);

   530 #endif

   531         }

   532     };

   533     template<class T>

   534     struct vector_norm_inf:

   535         public vector_scalar_real_unary_functor<T> {

   536         typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;

   537         typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;

   538         typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;

   539         typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;

   540         typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;

   542         template<class E>

   543         static BOOST_UBLAS_INLINE

   544         result_type apply (const vector_expression<E> &e) {

   545             real_type t = real_type ();

   546             size_type size (e ().size ());

   547             for (size_type i = 0; i < size; ++ i) {

   548                 real_type u (type_traits<value_type>::norm_inf (e () (i)));

   549                 if (u > t)

   550                     t = u;

   551             }

   552             return t;

   553         }

   554         // Dense case

   555         template<class I>

   556         static BOOST_UBLAS_INLINE

   557         result_type apply (difference_type size, I it) {

   558             real_type t = real_type ();

   559             while (-- size >= 0) {

   560                 real_type u (type_traits<value_type>::norm_inf (*it));

   561                 if (u > t)

   562                     t = u;

   563                 ++ it;

   564             }

   565             return t;

   566         }

   567         // Sparse case

   568         template<class I>

   569         static BOOST_UBLAS_INLINE

   570         result_type apply (I it, const I &it_end) {

   571             real_type t = real_type ();

   572             while (it != it_end) {

   573                 real_type u (type_traits<value_type>::norm_inf (*it));

   574                 if (u > t)

   575                     t = u;

   576                 ++ it;

   577             }

   578             return t;

   579         }

   580     };

   582     // Unary returning index

   583     template<class T>

   584     struct vector_scalar_index_unary_functor {

   585         typedef std::size_t size_type;

   586         typedef std::ptrdiff_t difference_type;

   587         typedef T value_type;

   588         typedef typename type_traits<T>::real_type real_type;

   589         typedef size_type result_type;

   590     };

   592     template<class T>

   593     struct vector_index_norm_inf:

   594         public vector_scalar_index_unary_functor<T> {

   595         typedef typename vector_scalar_index_unary_functor<T>::size_type size_type;

   596         typedef typename vector_scalar_index_unary_functor<T>::difference_type difference_type;

   597         typedef typename vector_scalar_index_unary_functor<T>::value_type value_type;

   598         typedef typename vector_scalar_index_unary_functor<T>::real_type real_type;

   599         typedef typename vector_scalar_index_unary_functor<T>::result_type result_type;

   601         template<class E>

   602         static BOOST_UBLAS_INLINE

   603         result_type apply (const vector_expression<E> &e) {

   604             // ISSUE For CBLAS compatibility return 0 index in empty case

   605             result_type i_norm_inf (0);

   606             real_type t = real_type ();

   607             size_type size (e ().size ());

   608             for (size_type i = 0; i < size; ++ i) {

   609                 real_type u (type_traits<value_type>::norm_inf (e () (i)));

   610                 if (u > t) {

   611                     i_norm_inf = i;

   612                     t = u;

   613                 }

   614             }

   615             return i_norm_inf;

   616         }

   617         // Dense case

   618         template<class I>

   619         static BOOST_UBLAS_INLINE

   620         result_type apply (difference_type size, I it) {

   621             // ISSUE For CBLAS compatibility return 0 index in empty case

   622             result_type i_norm_inf (0);

   623             real_type t = real_type ();

   624             while (-- size >= 0) {

   625                 real_type u (type_traits<value_type>::norm_inf (*it));

   626                 if (u > t) {

   627                     i_norm_inf = it.index ();

   628                     t = u;

   629                 }

   630                 ++ it;

   631             }

   632             return i_norm_inf;

   633         }

   634         // Sparse case

   635         template<class I>

   636         static BOOST_UBLAS_INLINE

   637         result_type apply (I it, const I &it_end) {

   638             // ISSUE For CBLAS compatibility return 0 index in empty case

   639             result_type i_norm_inf (0);

   640             real_type t = real_type ();

   641             while (it != it_end) {

   642                 real_type u (type_traits<value_type>::norm_inf (*it));

   643                 if (u > t) {

   644                     i_norm_inf = it.index ();

   645                     t = u;

   646                 }

   647                 ++ it;

   648             }

   649             return i_norm_inf;

   650         }

   651     };

   653     // Binary returning scalar

   654     template<class T1, class T2, class TR>

   655     struct vector_scalar_binary_functor {

   656         typedef std::size_t size_type;

   657         typedef std::ptrdiff_t difference_type;

   658         typedef TR value_type;

   659         typedef TR result_type;

   660     };

   662     template<class T1, class T2, class TR>

   663     struct vector_inner_prod:

   664         public vector_scalar_binary_functor<T1, T2, TR> {

   665         typedef typename vector_scalar_binary_functor<T1, T2, TR>::size_type size_type ;

   666         typedef typename vector_scalar_binary_functor<T1, T2, TR>::difference_type difference_type;

   667         typedef typename vector_scalar_binary_functor<T1, T2, TR>::value_type value_type;

   668         typedef typename vector_scalar_binary_functor<T1, T2, TR>::result_type result_type;

   670         template<class C1, class C2>

   671         static BOOST_UBLAS_INLINE

   672         result_type apply (const vector_container<C1> &c1,

   673                            const vector_container<C2> &c2) {

   674 #ifdef BOOST_UBLAS_USE_SIMD

   675             using namespace raw;

   676             size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));

   677             const T1 *data1 = data_const (c1 ());

   678             const T2 *data2 = data_const (c2 ());

   679             size_type s1 = stride (c1 ());

   680             size_type s2 = stride (c2 ());

   681             result_type t = result_type (0);

   682             if (s1 == 1 && s2 == 1) {

   683                 for (size_type i = 0; i < size; ++ i)

   684                     t += data1 [i] * data2 [i];

   685             } else if (s2 == 1) {

   686                 for (size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)

   687                     t += data1 [i1] * data2 [i];

   688             } else if (s1 == 1) {

   689                 for (size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)

   690                     t += data1 [i] * data2 [i2];

   691             } else {

   692                 for (size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)

   693                     t += data1 [i1] * data2 [i2];

   694             }

   695             return t;

   696 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)

   697             return boost::numeric::bindings::atlas::dot (c1 (), c2 ());

   698 #else

   699             return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2));

   700 #endif

   701         }

   702         template<class E1, class E2>

   703         static BOOST_UBLAS_INLINE

   704         result_type apply (const vector_expression<E1> &e1,

   705                            const vector_expression<E2> &e2) {

   706             size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));

   707             result_type t = result_type (0);

   708 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE

   709             for (size_type i = 0; i < size; ++ i)

   710                 t += e1 () (i) * e2 () (i);

   711 #else

   712             size_type i (0);

   713             DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));

   714 #endif

   715             return t;

   716         }

   717         // Dense case

   718         template<class I1, class I2>

   719         static BOOST_UBLAS_INLINE

   720         result_type apply (difference_type size, I1 it1, I2 it2) {

   721             result_type t = result_type (0);

   722 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE

   723             while (-- size >= 0)

   724                 t += *it1 * *it2, ++ it1, ++ it2;

   725 #else

   726             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));

   727 #endif

   728             return t;

   729         }

   730         // Packed case

   731         template<class I1, class I2>

   732         static BOOST_UBLAS_INLINE

   733         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {

   734             result_type t = result_type (0);

   735             difference_type it1_size (it1_end - it1);

   736             difference_type it2_size (it2_end - it2);

   737             difference_type diff (0);

   738             if (it1_size > 0 && it2_size > 0)

   739                 diff = it2.index () - it1.index ();

   740             if (diff != 0) {

   741                 difference_type size = (std::min) (diff, it1_size);

   742                 if (size > 0) {

   743                     it1 += size;

   744                     it1_size -= size;

   745                     diff -= size;

   746                 }

   747                 size = (std::min) (- diff, it2_size);

   748                 if (size > 0) {

   749                     it2 += size;

   750                     it2_size -= size;

   751                     diff += size;

   752                 }

   753             }

   754             difference_type size ((std::min) (it1_size, it2_size));

   755             while (-- size >= 0)

   756                 t += *it1 * *it2, ++ it1, ++ it2;

   757             return t;

   758         }

   759         // Sparse case

   760         template<class I1, class I2>

   761         static BOOST_UBLAS_INLINE

   762         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {

   763             result_type t = result_type (0);

   764             if (it1 != it1_end && it2 != it2_end) {

   765                 size_type it1_index = it1.index (), it2_index = it2.index ();

   766                 while (true) {

   767                     difference_type compare = it1_index - it2_index;

   768                     if (compare == 0) {

   769                         t += *it1 * *it2, ++ it1, ++ it2;

   770                         if (it1 != it1_end && it2 != it2_end) {

   771                             it1_index = it1.index ();

   772                             it2_index = it2.index ();

   773                         } else

   774                             break;

   775                     } else if (compare < 0) {

   776                         increment (it1, it1_end, - compare);

   777                         if (it1 != it1_end)

   778                             it1_index = it1.index ();

   779                         else

   780                             break;

   781                     } else if (compare > 0) {

   782                         increment (it2, it2_end, compare);

   783                         if (it2 != it2_end)

   784                             it2_index = it2.index ();

   785                         else

   786                             break;

   787                     }

   788                 }

   789             }

   790             return t;

   791         }

   792     };

   794     // Matrix functors

   796     // Binary returning vector

   797     template<class T1, class T2, class TR>

   798     struct matrix_vector_binary_functor {

   799         typedef std::size_t size_type;

   800         typedef std::ptrdiff_t difference_type;

   801         typedef TR value_type;

   802         typedef TR result_type;

   803     };

   805     template<class T1, class T2, class TR>

   806     struct matrix_vector_prod1:

   807         public matrix_vector_binary_functor<T1, T2, TR> {

   808         typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;

   809         typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;

   810         typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;

   811         typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;

   813         template<class C1, class C2>

   814         static BOOST_UBLAS_INLINE

   815         result_type apply (const matrix_container<C1> &c1,

   816                            const vector_container<C2> &c2,

   817                            size_type i) {

   818 #ifdef BOOST_UBLAS_USE_SIMD

   819             using namespace raw;

   820             size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ());

   821             const T1 *data1 = data_const (c1 ()) + i * stride1 (c1 ());

   822             const T2 *data2 = data_const (c2 ());

   823             size_type s1 = stride2 (c1 ());

   824             size_type s2 = stride (c2 ());

   825             result_type t = result_type (0);

   826             if (s1 == 1 && s2 == 1) {

   827                 for (size_type j = 0; j < size; ++ j)

   828                     t += data1 [j] * data2 [j];

   829             } else if (s2 == 1) {

   830                 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)

   831                     t += data1 [j1] * data2 [j];

   832             } else if (s1 == 1) {

   833                 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)

   834                     t += data1 [j] * data2 [j2];

   835             } else {

   836                 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)

   837                     t += data1 [j1] * data2 [j2];

   838             }

   839             return t;

   840 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)

   841             return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ());

   842 #else

   843             return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i));

   844 #endif

   845         }

   846         template<class E1, class E2>

   847         static BOOST_UBLAS_INLINE

   848         result_type apply (const matrix_expression<E1> &e1,

   849                                  const vector_expression<E2> &e2,

   850                            size_type i) {

   851             size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());

   852             result_type t = result_type (0);

   853 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE

   854             for (size_type j = 0; j < size; ++ j)

   855                 t += e1 () (i, j) * e2 () (j);

   856 #else

   857             size_type j (0);

   858             DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j));

   859 #endif

   860             return t;

   861         }

   862         // Dense case

   863         template<class I1, class I2>

   864         static BOOST_UBLAS_INLINE

   865         result_type apply (difference_type size, I1 it1, I2 it2) {

   866             result_type t = result_type (0);

   867 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE

   868             while (-- size >= 0)

   869                 t += *it1 * *it2, ++ it1, ++ it2;

   870 #else

   871             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));

   872 #endif

   873             return t;

   874         }

   875         // Packed case

   876         template<class I1, class I2>

   877         static BOOST_UBLAS_INLINE

   878         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {

   879             result_type t = result_type (0);

   880             difference_type it1_size (it1_end - it1);

   881             difference_type it2_size (it2_end - it2);

   882             difference_type diff (0);

   883             if (it1_size > 0 && it2_size > 0)

   884                 diff = it2.index () - it1.index2 ();

   885             if (diff != 0) {

   886                 difference_type size = (std::min) (diff, it1_size);

   887                 if (size > 0) {

   888                     it1 += size;

   889                     it1_size -= size;

   890                     diff -= size;

   891                 }

   892                 size = (std::min) (- diff, it2_size);

   893                 if (size > 0) {

   894                     it2 += size;

   895                     it2_size -= size;

   896                     diff += size;

   897                 }

   898             }

   899             difference_type size ((std::min) (it1_size, it2_size));

   900             while (-- size >= 0)

   901                 t += *it1 * *it2, ++ it1, ++ it2;

   902             return t;

   903         }

   904         // Sparse case

   905         template<class I1, class I2>

   906         static BOOST_UBLAS_INLINE

   907         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,

   908                                  sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {

   909             result_type t = result_type (0);

   910             if (it1 != it1_end && it2 != it2_end) {

   911                 size_type it1_index = it1.index2 (), it2_index = it2.index ();

   912                 while (true) {

   913                     difference_type compare = it1_index - it2_index;

   914                     if (compare == 0) {

   915                         t += *it1 * *it2, ++ it1, ++ it2;

   916                         if (it1 != it1_end && it2 != it2_end) {

   917                             it1_index = it1.index2 ();

   918                             it2_index = it2.index ();

   919                         } else

   920                             break;

   921                     } else if (compare < 0) {

   922                         increment (it1, it1_end, - compare);

   923                         if (it1 != it1_end)

   924                             it1_index = it1.index2 ();

   925                         else

   926                             break;

   927                     } else if (compare > 0) {

   928                         increment (it2, it2_end, compare);

   929                         if (it2 != it2_end)

   930                             it2_index = it2.index ();

   931                         else

   932                             break;

   933                     }

   934                 }

   935             }

   936             return t;

   937         }

   938         // Sparse packed case

   939         template<class I1, class I2>

   940         static BOOST_UBLAS_INLINE

   941         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,

   942                                  sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {

   943             result_type t = result_type (0);

   944             while (it1 != it1_end) {

   945                 t += *it1 * it2 () (it1.index2 ());

   946                 ++ it1;

   947             }

   948             return t;

   949         }

   950         // Packed sparse case

   951         template<class I1, class I2>

   952         static BOOST_UBLAS_INLINE

   953         result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,

   954                                  packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {

   955             result_type t = result_type (0);

   956             while (it2 != it2_end) {

   957                 t += it1 () (it1.index1 (), it2.index ()) * *it2;

   958                 ++ it2;

   959             }

   960             return t;

   961         }

   962         // Another dispatcher

   963         template<class I1, class I2>

   964         static BOOST_UBLAS_INLINE

   965         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,

   966                                  sparse_bidirectional_iterator_tag) {

   967             typedef typename I1::iterator_category iterator1_category;

   968             typedef typename I2::iterator_category iterator2_category;

   969             return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());

   970         }

   971     };

   973     template<class T1, class T2, class TR>

   974     struct matrix_vector_prod2:

   975         public matrix_vector_binary_functor<T1, T2, TR> {

   976         typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;

   977         typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;

   978         typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;

   979         typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;

   981         template<class C1, class C2>

   982         static BOOST_UBLAS_INLINE

   983         result_type apply (const vector_container<C1> &c1,

   984                            const matrix_container<C2> &c2,

   985                            size_type i) {

   986 #ifdef BOOST_UBLAS_USE_SIMD

   987             using namespace raw;

   988             size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ());

   989             const T1 *data1 = data_const (c1 ());

   990             const T2 *data2 = data_const (c2 ()) + i * stride2 (c2 ());

   991             size_type s1 = stride (c1 ());

   992             size_type s2 = stride1 (c2 ());

   993             result_type t = result_type (0);

   994             if (s1 == 1 && s2 == 1) {

   995                 for (size_type j = 0; j < size; ++ j)

   996                     t += data1 [j] * data2 [j];

   997             } else if (s2 == 1) {

   998                 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)

   999                     t += data1 [j1] * data2 [j];

  1000             } else if (s1 == 1) {

  1001                 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)

  1002                     t += data1 [j] * data2 [j2];

  1003             } else {

  1004                 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)

  1005                     t += data1 [j1] * data2 [j2];

  1006             }

  1007             return t;

  1008 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)

  1009             return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i));

  1010 #else

  1011             return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));

  1012 #endif

  1013         }

  1014         template<class E1, class E2>

  1015         static BOOST_UBLAS_INLINE

  1016         result_type apply (const vector_expression<E1> &e1,

  1017                                  const matrix_expression<E2> &e2,

  1018                            size_type i) {

  1019             size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());

  1020             result_type t = result_type (0);

  1021 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE

  1022             for (size_type j = 0; j < size; ++ j)

  1023                 t += e1 () (j) * e2 () (j, i);

  1024 #else

  1025             size_type j (0);

  1026             DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j));

  1027 #endif

  1028             return t;

  1029         }

  1030         // Dense case

  1031         template<class I1, class I2>

  1032         static BOOST_UBLAS_INLINE

  1033         result_type apply (difference_type size, I1 it1, I2 it2) {

  1034             result_type t = result_type (0);

  1035 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE

  1036             while (-- size >= 0)

  1037                 t += *it1 * *it2, ++ it1, ++ it2;

  1038 #else

  1039             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));

  1040 #endif

  1041             return t;

  1042         }

  1043         // Packed case

  1044         template<class I1, class I2>

  1045         static BOOST_UBLAS_INLINE

  1046         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {

  1047             result_type t = result_type (0);

  1048             difference_type it1_size (it1_end - it1);

  1049             difference_type it2_size (it2_end - it2);

  1050             difference_type diff (0);

  1051             if (it1_size > 0 && it2_size > 0)

  1052                 diff = it2.index1 () - it1.index ();

  1053             if (diff != 0) {

  1054                 difference_type size = (std::min) (diff, it1_size);

  1055                 if (size > 0) {

  1056                     it1 += size;

  1057                     it1_size -= size;

  1058                     diff -= size;

  1059                 }

  1060                 size = (std::min) (- diff, it2_size);

  1061                 if (size > 0) {

  1062                     it2 += size;

  1063                     it2_size -= size;

  1064                     diff += size;

  1065                 }

  1066             }

  1067             difference_type size ((std::min) (it1_size, it2_size));

  1068             while (-- size >= 0)

  1069                 t += *it1 * *it2, ++ it1, ++ it2;

  1070             return t;

  1071         }

  1072         // Sparse case

  1073         template<class I1, class I2>

  1074         static BOOST_UBLAS_INLINE

  1075         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,

  1076                                  sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {

  1077             result_type t = result_type (0);

  1078             if (it1 != it1_end && it2 != it2_end) {

  1079                 size_type it1_index = it1.index (), it2_index = it2.index1 ();

  1080                 while (true) {

  1081                     difference_type compare = it1_index - it2_index;

  1082                     if (compare == 0) {

  1083                         t += *it1 * *it2, ++ it1, ++ it2;

  1084                         if (it1 != it1_end && it2 != it2_end) {

  1085                             it1_index = it1.index ();

  1086                             it2_index = it2.index1 ();

  1087                         } else

  1088                             break;

  1089                     } else if (compare < 0) {

  1090                         increment (it1, it1_end, - compare);

  1091                         if (it1 != it1_end)

  1092                             it1_index = it1.index ();

  1093                         else

  1094                             break;

  1095                     } else if (compare > 0) {

  1096                         increment (it2, it2_end, compare);

  1097                         if (it2 != it2_end)

  1098                             it2_index = it2.index1 ();

  1099                         else

  1100                             break;

  1101                     }

  1102                 }

  1103             }

  1104             return t;

  1105         }

  1106         // Packed sparse case

  1107         template<class I1, class I2>

  1108         static BOOST_UBLAS_INLINE

  1109         result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,

  1110                                  packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {

  1111             result_type t = result_type (0);

  1112             while (it2 != it2_end) {

  1113                 t += it1 () (it2.index1 ()) * *it2;

  1114                 ++ it2;

  1115             }

  1116             return t;

  1117         }

  1118         // Sparse packed case

  1119         template<class I1, class I2>

  1120         static BOOST_UBLAS_INLINE

  1121         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,

  1122                                  sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {

  1123             result_type t = result_type (0);

  1124             while (it1 != it1_end) {

  1125                 t += *it1 * it2 () (it1.index (), it2.index2 ());

  1126                 ++ it1;

  1127             }

  1128             return t;

  1129         }

  1130         // Another dispatcher

  1131         template<class I1, class I2>

  1132         static BOOST_UBLAS_INLINE

  1133         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,

  1134                                  sparse_bidirectional_iterator_tag) {

  1135             typedef typename I1::iterator_category iterator1_category;

  1136             typedef typename I2::iterator_category iterator2_category;

  1137             return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());

  1138         }

  1139     };

  1141     // Binary returning matrix

  1142     template<class T1, class T2, class TR>

  1143     struct matrix_matrix_binary_functor {

  1144         typedef std::size_t size_type;

  1145         typedef std::ptrdiff_t difference_type;

  1146         typedef TR value_type;

  1147         typedef TR result_type;

  1148     };

  1150     template<class T1, class T2, class TR>

  1151     struct matrix_matrix_prod:

  1152         public matrix_matrix_binary_functor<T1, T2, TR> {

  1153         typedef typename matrix_matrix_binary_functor<T1, T2, TR>::size_type size_type;

  1154         typedef typename matrix_matrix_binary_functor<T1, T2, TR>::difference_type difference_type;

  1155         typedef typename matrix_matrix_binary_functor<T1, T2, TR>::value_type value_type;

  1156         typedef typename matrix_matrix_binary_functor<T1, T2, TR>::result_type result_type;

  1158         template<class C1, class C2>

  1159         static BOOST_UBLAS_INLINE

  1160         result_type apply (const matrix_container<C1> &c1,

  1161                            const matrix_container<C2> &c2,

  1162                            size_type i, size_type j) {

  1163 #ifdef BOOST_UBLAS_USE_SIMD

  1164             using namespace raw;

  1165             size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ());

  1166             const T1 *data1 = data_const (c1 ()) + i * stride1 (c1 ());

  1167             const T2 *data2 = data_const (c2 ()) + j * stride2 (c2 ());

  1168             size_type s1 = stride2 (c1 ());

  1169             size_type s2 = stride1 (c2 ());

  1170             result_type t = result_type (0);

  1171             if (s1 == 1 && s2 == 1) {

  1172                 for (size_type k = 0; k < size; ++ k)

  1173                     t += data1 [k] * data2 [k];

  1174             } else if (s2 == 1) {

  1175                 for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1)

  1176                     t += data1 [k1] * data2 [k];

  1177             } else if (s1 == 1) {

  1178                 for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2)

  1179                     t += data1 [k] * data2 [k2];

  1180             } else {

  1181                 for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2)

  1182                     t += data1 [k1] * data2 [k2];

  1183             }

  1184             return t;

  1185 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)

  1186             return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j));

  1187 #else

  1188             return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));

  1189 #endif

  1190         }

  1191         template<class E1, class E2>

  1192         static BOOST_UBLAS_INLINE

  1193         result_type apply (const matrix_expression<E1> &e1,

  1194                                  const matrix_expression<E2> &e2,

  1195                            size_type i, size_type j) {

  1196             size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());

  1197             result_type t = result_type (0);

  1198 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE

  1199             for (size_type k = 0; k < size; ++ k)

  1200                 t += e1 () (i, k) * e2 () (k, j);

  1201 #else

  1202             size_type k (0);

  1203             DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k));

  1204 #endif

  1205             return t;

  1206         }

  1207         // Dense case

  1208         template<class I1, class I2>

  1209         static BOOST_UBLAS_INLINE

  1210         result_type apply (difference_type size, I1 it1, I2 it2) {

  1211             result_type t = result_type (0);

  1212 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE

  1213             while (-- size >= 0)

  1214                 t += *it1 * *it2, ++ it1, ++ it2;

  1215 #else

  1216             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));

  1217 #endif

  1218             return t;

  1219         }

  1220         // Packed case

  1221         template<class I1, class I2>

  1222         static BOOST_UBLAS_INLINE

  1223         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) {

  1224             result_type t = result_type (0);

  1225             difference_type it1_size (it1_end - it1);

  1226             difference_type it2_size (it2_end - it2);

  1227             difference_type diff (0);

  1228             if (it1_size > 0 && it2_size > 0)

  1229                 diff = it2.index1 () - it1.index2 ();

  1230             if (diff != 0) {

  1231                 difference_type size = (std::min) (diff, it1_size);

  1232                 if (size > 0) {

  1233                     it1 += size;

  1234                     it1_size -= size;

  1235                     diff -= size;

  1236                 }

  1237                 size = (std::min) (- diff, it2_size);

  1238                 if (size > 0) {

  1239                     it2 += size;

  1240                     it2_size -= size;

  1241                     diff += size;

  1242                 }

  1243             }

  1244             difference_type size ((std::min) (it1_size, it2_size));

  1245             while (-- size >= 0)

  1246                 t += *it1 * *it2, ++ it1, ++ it2;

  1247             return t;

  1248         }

  1249         // Sparse case

  1250         template<class I1, class I2>

  1251         static BOOST_UBLAS_INLINE

  1252         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {

  1253             result_type t = result_type (0);

  1254             if (it1 != it1_end && it2 != it2_end) {

  1255                 size_type it1_index = it1.index2 (), it2_index = it2.index1 ();

  1256                 while (true) {

  1257                     difference_type compare = it1_index - it2_index;

  1258                     if (compare == 0) {

  1259                         t += *it1 * *it2, ++ it1, ++ it2;

  1260                         if (it1 != it1_end && it2 != it2_end) {

  1261                             it1_index = it1.index2 ();

  1262                             it2_index = it2.index1 ();

  1263                         } else

  1264                             break;

  1265                     } else if (compare < 0) {

  1266                         increment (it1, it1_end, - compare);

  1267                         if (it1 != it1_end)

  1268                             it1_index = it1.index2 ();

  1269                         else

  1270                             break;

  1271                     } else if (compare > 0) {

  1272                         increment (it2, it2_end, compare);

  1273                         if (it2 != it2_end)

  1274                             it2_index = it2.index1 ();

  1275                         else

  1276                             break;

  1277                     }

  1278                 }

  1279             }

  1280             return t;

  1281         }

  1282     };

  1284     // Unary returning scalar norm

  1285     template<class T>

  1286     struct matrix_scalar_real_unary_functor {

  1287         typedef std::size_t size_type;

  1288         typedef std::ptrdiff_t difference_type;

  1289         typedef T value_type;

  1290         typedef typename type_traits<T>::real_type real_type;

  1291         typedef real_type result_type;

  1292     };

  1294     template<class T>

  1295     struct matrix_norm_1:

  1296         public matrix_scalar_real_unary_functor<T> {

  1297         typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;

  1298         typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;

  1299         typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;

  1300         typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;

  1301         typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;

  1303         template<class E>

  1304         static BOOST_UBLAS_INLINE

  1305         result_type apply (const matrix_expression<E> &e) {

  1306             real_type t = real_type ();

  1307             size_type size2 (e ().size2 ());

  1308             for (size_type j = 0; j < size2; ++ j) {

  1309                 real_type u = real_type ();

  1310                 size_type size1 (e ().size1 ());

  1311                 for (size_type i = 0; i < size1; ++ i) {

  1312                     real_type v (type_traits<value_type>::norm_1 (e () (i, j)));

  1313                     u += v;

  1314                 }

  1315                 if (u > t)

  1316                     t = u;

  1317             }

  1318             return t;

  1319         }

  1320     };

  1321     template<class T>

  1322     struct matrix_norm_frobenius:

  1323         public matrix_scalar_real_unary_functor<T> {

  1324         typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;

  1325         typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;

  1326         typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;

  1327         typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;

  1328         typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;

  1330         template<class E>

  1331         static BOOST_UBLAS_INLINE

  1332         result_type apply (const matrix_expression<E> &e) {

  1333             real_type t = real_type ();

  1334             size_type size1 (e ().size1 ());

  1335             for (size_type i = 0; i < size1; ++ i) {

  1336                 size_type size2 (e ().size2 ());

  1337                 for (size_type j = 0; j < size2; ++ j) {

  1338                     real_type u (type_traits<value_type>::norm_2 (e () (i, j)));

  1339                     t +=  u * u;

  1340                 }

  1341             }

  1342             return type_traits<real_type>::type_sqrt (t);

  1343         }

  1344     };

  1345     template<class T>

  1346     struct matrix_norm_inf:

  1347         public matrix_scalar_real_unary_functor<T> {

  1348         typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;

  1349         typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;

  1350         typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;

  1351         typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;

  1352         typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;

  1354         template<class E>

  1355         static BOOST_UBLAS_INLINE

  1356         result_type apply (const matrix_expression<E> &e) {

  1357             real_type t = real_type ();

  1358             size_type size1 (e ().size1 ());

  1359             for (size_type i = 0; i < size1; ++ i) {

  1360                 real_type u = real_type ();

  1361                 size_type size2 (e ().size2 ());

  1362                 for (size_type j = 0; j < size2; ++ j) {

  1363                     real_type v (type_traits<value_type>::norm_inf (e () (i, j)));

  1364                     u += v;

  1365                 }

  1366                 if (u > t)

  1367                     t = u;

  1368             }

  1369             return t;

  1370         }

  1371     };

  1373     // This functor computes the address translation

  1374     // matrix [i] [j] -> storage [i * size2 + j]

  1375     template <class Z, class D>

  1376     struct basic_row_major {

  1377         typedef Z size_type;

  1378         typedef D difference_type;

  1379         typedef row_major_tag orientation_category;

  1381         static

  1382         BOOST_UBLAS_INLINE

  1383         size_type storage_size (size_type size1, size_type size2) {

  1384             // Guard against size_type overflow

  1385             BOOST_UBLAS_CHECK (size2 == 0 || size1 <= (std::numeric_limits<size_type>::max) () / size2, bad_size ());

  1386             return size1 * size2;

  1387         }

  1389         // Indexing

  1390         static

  1391         BOOST_UBLAS_INLINE

  1392         size_type element (size_type i, size_type size1, size_type j, size_type size2) {

  1393             BOOST_UBLAS_CHECK (i < size1, bad_index ());

  1394             BOOST_UBLAS_CHECK (j < size2, bad_index ());

  1395             detail::ignore_unused_variable_warning(size1);

  1396             // Guard against size_type overflow

  1397             BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size2, bad_index ());

  1398             return i * size2 + j;

  1399         }

  1400         static

  1401         BOOST_UBLAS_INLINE

  1402         size_type address (size_type i, size_type size1, size_type j, size_type size2) {

  1403             BOOST_UBLAS_CHECK (i <= size1, bad_index ());

  1404             BOOST_UBLAS_CHECK (j <= size2, bad_index ());

  1405             // Guard against size_type overflow - address may be size2 past end of storage

  1406             BOOST_UBLAS_CHECK (size2 == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size2, bad_index ());

  1407             detail::ignore_unused_variable_warning(size1);

  1408             return i * size2 + j;

  1409         }

  1410         static

  1411         BOOST_UBLAS_INLINE

  1412         difference_type distance1 (difference_type k, size_type /* size1 */, size_type size2) {

  1413             return size2 != 0 ? k / size2 : 0;

  1414         }

  1415         static

  1416         BOOST_UBLAS_INLINE

  1417         difference_type distance2 (difference_type k, size_type /* size1 */, size_type /* size2 */) {

  1418             return k;

  1419         }

  1420         static

  1421         BOOST_UBLAS_INLINE

  1422         size_type index1 (difference_type k, size_type /* size1 */, size_type size2) {

  1423             return size2 != 0 ? k / size2 : 0;

  1424         }

  1425         static

  1426         BOOST_UBLAS_INLINE

  1427         size_type index2 (difference_type k, size_type /* size1 */, size_type size2) {

  1428             return size2 != 0 ? k % size2 : 0;

  1429         }

  1430         static

  1431         BOOST_UBLAS_INLINE

  1432         bool fast1 () {

  1433             return false;

  1434         }

  1435         static

  1436         BOOST_UBLAS_INLINE

  1437         size_type one1 (size_type /* size1 */, size_type size2) {

  1438             return size2;

  1439         }

  1440         static

  1441         BOOST_UBLAS_INLINE

  1442         bool fast2 () {

  1443             return true;

  1444         }

  1445         static

  1446         BOOST_UBLAS_INLINE

  1447         size_type one2 (size_type /* size1 */, size_type /* size2 */) {

  1448             return 1;

  1449         }

  1451         static

  1452         BOOST_UBLAS_INLINE

  1453         size_type triangular_size (size_type size1, size_type size2) {

  1454             size_type size = (std::max) (size1, size2);

  1455             // Guard against size_type overflow - simplified

  1456             BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());

  1457             return ((size + 1) * size) / 2;

  1458         }

  1459         static

  1460         BOOST_UBLAS_INLINE

  1461         size_type lower_element (size_type i, size_type size1, size_type j, size_type size2) {

  1462             BOOST_UBLAS_CHECK (i < size1, bad_index ());

  1463             BOOST_UBLAS_CHECK (j < size2, bad_index ());

  1464             BOOST_UBLAS_CHECK (i >= j, bad_index ());

  1465             detail::ignore_unused_variable_warning(size1);

  1466             detail::ignore_unused_variable_warning(size2);

  1467             // FIXME size_type overflow

  1468             // sigma_i (i + 1) = (i + 1) * i / 2

  1469             // i = 0 1 2 3, sigma = 0 1 3 6

  1470             return ((i + 1) * i) / 2 + j;

  1471         }

  1472         static

  1473         BOOST_UBLAS_INLINE

  1474         size_type upper_element (size_type i, size_type size1, size_type j, size_type size2) {

  1475             BOOST_UBLAS_CHECK (i < size1, bad_index ());

  1476             BOOST_UBLAS_CHECK (j < size2, bad_index ());

  1477             BOOST_UBLAS_CHECK (i <= j, bad_index ());

  1478             // FIXME size_type overflow

  1479             // sigma_i (size - i) = size * i - i * (i - 1) / 2

  1480             // i = 0 1 2 3, sigma = 0 4 7 9

  1481             return (i * (2 * (std::max) (size1, size2) - i + 1)) / 2 + j - i;

  1482         }

  1484         static

  1485         BOOST_UBLAS_INLINE

  1486         size_type element1 (size_type i, size_type size1, size_type /* j */, size_type /* size2 */) {

  1487             BOOST_UBLAS_CHECK (i < size1, bad_index ());

  1488             detail::ignore_unused_variable_warning(size1);

  1489             return i;

  1490         }

  1491         static

  1492         BOOST_UBLAS_INLINE

  1493         size_type element2 (size_type /* i */, size_type /* size1 */, size_type j, size_type size2) {

  1494             BOOST_UBLAS_CHECK (j < size2, bad_index ());

  1495             detail::ignore_unused_variable_warning(size2);

  1496             return j;

  1497         }

  1498         static

  1499         BOOST_UBLAS_INLINE

  1500         size_type address1 (size_type i, size_type size1, size_type /* j */, size_type /* size2 */) {

  1501             BOOST_UBLAS_CHECK (i <= size1, bad_index ());

  1502             detail::ignore_unused_variable_warning(size1);

  1503             return i;

  1504         }

  1505         static

  1506         BOOST_UBLAS_INLINE

  1507         size_type address2 (size_type /* i */, size_type /* size1 */, size_type j, size_type size2) {

  1508             BOOST_UBLAS_CHECK (j <= size2, bad_index ());

  1509             detail::ignore_unused_variable_warning(size2);

  1510             return j;

  1511         }

  1512         static

  1513         BOOST_UBLAS_INLINE

  1514         size_type index1 (size_type index1, size_type /* index2 */) {

  1515             return index1;

  1516         }

  1517         static

  1518         BOOST_UBLAS_INLINE

  1519         size_type index2 (size_type /* index1 */, size_type index2) {

  1520             return index2;

  1521         }

  1522         static

  1523         BOOST_UBLAS_INLINE

  1524         size_type size1 (size_type size1, size_type /* size2 */) {

  1525             return size1;

  1526         }

  1527         static

  1528         BOOST_UBLAS_INLINE

  1529         size_type size2 (size_type /* size1 */, size_type size2) {

  1530             return size2;

  1531         }

  1533         // Iterating

  1534         template<class I>

  1535         static

  1536         BOOST_UBLAS_INLINE

  1537         void increment1 (I &it, size_type /* size1 */, size_type size2) {

  1538             it += size2;

  1539         }

  1540         template<class I>

  1541         static

  1542         BOOST_UBLAS_INLINE

  1543         void decrement1 (I &it, size_type /* size1 */, size_type size2) {

  1544             it -= size2;

  1545         }

  1546         template<class I>

  1547         static

  1548         BOOST_UBLAS_INLINE

  1549         void increment2 (I &it, size_type /* size1 */, size_type /* size2 */) {

  1550             ++ it;

  1551         }

  1552         template<class I>

  1553         static

  1554         BOOST_UBLAS_INLINE

  1555         void decrement2 (I &it, size_type /* size1 */, size_type /* size2 */) {

  1556             -- it;

  1557         }

  1558     };

  1560     // This functor computes the address translation

  1561     // matrix [i] [j] -> storage [i + j * size1]

  1562     template <class Z, class D>

  1563     struct basic_column_major {

  1564         typedef Z size_type;

  1565         typedef D difference_type;

  1566         typedef column_major_tag orientation_category;

  1568         static

  1569         BOOST_UBLAS_INLINE

  1570         size_type storage_size (size_type size1, size_type size2) {

  1571             // Guard against size_type overflow

  1572             BOOST_UBLAS_CHECK (size1 == 0 || size2 <= (std::numeric_limits<size_type>::max) () / size1, bad_size ());

  1573             return size1 * size2;

  1574         }

  1576         // Indexing

  1577         static

  1578         BOOST_UBLAS_INLINE

  1579         size_type element (size_type i, size_type size1, size_type j, size_type size2) {

  1580             BOOST_UBLAS_CHECK (i < size1, bad_index ());

  1581             BOOST_UBLAS_CHECK (j < size2, bad_index ());

  1582             detail::ignore_unused_variable_warning(size2);

  1583             // Guard against size_type overflow

  1584             BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size1, bad_index ());

  1585             return i + j * size1;

  1586         }

  1587         static

  1588         BOOST_UBLAS_INLINE

  1589         size_type address (size_type i, size_type size1, size_type j, size_type size2) {

  1590             BOOST_UBLAS_CHECK (i <= size1, bad_index ());

  1591             BOOST_UBLAS_CHECK (j <= size2, bad_index ());

  1592             detail::ignore_unused_variable_warning(size2);

  1593             // Guard against size_type overflow - address may be size1 past end of storage

  1594             BOOST_UBLAS_CHECK (size1 == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size1, bad_index ());

  1595             return i + j * size1;

  1596         }

  1597         static

  1598         BOOST_UBLAS_INLINE

  1599         difference_type distance1 (difference_type k, size_type /* size1 */, size_type /* size2 */) {

  1600             return k;

  1601         }

  1602         static

  1603         BOOST_UBLAS_INLINE

  1604         difference_type distance2 (difference_type k, size_type size1, size_type /* size2 */) {

  1605             return size1 != 0 ? k / size1 : 0;

  1606         }

  1607         static

  1608         BOOST_UBLAS_INLINE

  1609         size_type index1 (difference_type k, size_type size1, size_type /* size2 */) {

  1610             return size1 != 0 ? k % size1 : 0;

  1611         }

  1612         static

  1613         BOOST_UBLAS_INLINE

  1614         size_type index2 (difference_type k, size_type size1, size_type /* size2 */) {

  1615             return size1 != 0 ? k / size1 : 0;

  1616         }

  1617         static

  1618         BOOST_UBLAS_INLINE

  1619         bool fast1 () {

  1620             return true;

  1621         }

  1622         static

  1623         BOOST_UBLAS_INLINE

  1624         size_type one1 (size_type /* size1 */, size_type /* size2 */) {

  1625             return 1;

  1626         }

  1627         static

  1628         BOOST_UBLAS_INLINE

  1629         bool fast2 () {

  1630             return false;

  1631         }

  1632         static

  1633         BOOST_UBLAS_INLINE

  1634         size_type one2 (size_type size1, size_type /* size2 */) {

  1635             return size1;

  1636         }

  1638         static

  1639         BOOST_UBLAS_INLINE

  1640         size_type triangular_size (size_type size1, size_type size2) {

  1641             size_type size = (std::max) (size1, size2);

  1642             // Guard against size_type overflow - simplified

  1643             BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());

  1644             return ((size + 1) * size) / 2;

  1645         }

  1646         static

  1647         BOOST_UBLAS_INLINE

  1648         size_type lower_element (size_type i, size_type size1, size_type j, size_type size2) {

  1649             BOOST_UBLAS_CHECK (i < size1, bad_index ());

  1650             BOOST_UBLAS_CHECK (j < size2, bad_index ());

  1651             BOOST_UBLAS_CHECK (i >= j, bad_index ());

  1652             // FIXME size_type overflow

  1653             // sigma_j (size - j) = size * j - j * (j - 1) / 2

  1654             // j = 0 1 2 3, sigma = 0 4 7 9

  1655             return i - j + (j * (2 * (std::max) (size1, size2) - j + 1)) / 2;

  1656         }

  1657         static

  1658         BOOST_UBLAS_INLINE

  1659         size_type upper_element (size_type i, size_type size1, size_type j, size_type size2) {

  1660             BOOST_UBLAS_CHECK (i < size1, bad_index ());

  1661             BOOST_UBLAS_CHECK (j < size2, bad_index ());

  1662             BOOST_UBLAS_CHECK (i <= j, bad_index ());

  1663             // FIXME size_type overflow

  1664             // sigma_j (j + 1) = (j + 1) * j / 2

  1665             // j = 0 1 2 3, sigma = 0 1 3 6

  1666             return i + ((j + 1) * j) / 2;

  1667         }

  1669         static

  1670         BOOST_UBLAS_INLINE

  1671         size_type element1 (size_type /* i */, size_type /* size1 */, size_type j, size_type size2) {

  1672             BOOST_UBLAS_CHECK (j < size2, bad_index ());

  1673             detail::ignore_unused_variable_warning(size2);

  1674             return j;

  1675         }

  1676         static

  1677         BOOST_UBLAS_INLINE

  1678         size_type element2 (size_type i, size_type size1, size_type /* j */, size_type /* size2 */) {

  1679             BOOST_UBLAS_CHECK (i < size1, bad_index ());

  1680             detail::ignore_unused_variable_warning(size1);

  1681             return i;

  1682         }

  1683         static

  1684         BOOST_UBLAS_INLINE

  1685         size_type address1 (size_type /* i */, size_type /* size1 */, size_type j, size_type size2) {

  1686             BOOST_UBLAS_CHECK (j <= size2, bad_index ());

  1687             detail::ignore_unused_variable_warning(size2);

  1688             return j;

  1689         }

  1690         static

  1691         BOOST_UBLAS_INLINE

  1692         size_type address2 (size_type i, size_type size1, size_type /* j */, size_type /* size2 */) {

  1693             BOOST_UBLAS_CHECK (i <= size1, bad_index ());

  1694             detail::ignore_unused_variable_warning(size1);

  1695             return i;

  1696         }

  1697         static

  1698         BOOST_UBLAS_INLINE

  1699         size_type index1 (size_type /* index1 */, size_type index2) {

  1700             return index2;

  1701         }

  1702         static

  1703         BOOST_UBLAS_INLINE

  1704         size_type index2 (size_type index1, size_type /* index2 */) {

  1705             return index1;

  1706         }

  1707         static

  1708         BOOST_UBLAS_INLINE

  1709         size_type size1 (size_type /* size1 */, size_type size2) {

  1710             return size2;

  1711         }

  1712         static

  1713         BOOST_UBLAS_INLINE

  1714         size_type size2 (size_type size1, size_type /* size2 */) {

  1715             return size1;

  1716         }

  1718         // Iterating

  1719         template<class I>

  1720         static

  1721         BOOST_UBLAS_INLINE

  1722         void increment1 (I &it, size_type /* size1 */, size_type /* size2 */) {

  1723             ++ it;

  1724         }

  1725         template<class I>

  1726         static

  1727         BOOST_UBLAS_INLINE

  1728         void decrement1 (I &it, size_type /* size1 */, size_type /* size2 */) {

  1729             -- it;

  1730         }

  1731         template<class I>

  1732         static

  1733         BOOST_UBLAS_INLINE

  1734         void increment2 (I &it, size_type size1, size_type /* size2 */) {

  1735             it += size1;

  1736         }

  1737         template<class I>

  1738         static

  1739         BOOST_UBLAS_INLINE

  1740         void decrement2 (I &it, size_type size1, size_type /* size2 */) {

  1741             it -= size1;

  1742         }

  1743     };

  1746     template <class Z>

  1747     struct basic_full {

  1748         typedef Z size_type;

  1750         template<class L>

  1751         static

  1752         BOOST_UBLAS_INLINE

  1753         size_type packed_size (size_type size1, size_type size2) {

  1754             return L::storage_size (size1, size2);

  1755         }

  1757         static

  1758         BOOST_UBLAS_INLINE

  1759         bool zero (size_type /* i */, size_type /* j */) {

  1760             return false;

  1761         }

  1762         static

  1763         BOOST_UBLAS_INLINE

  1764         bool one (size_type /* i */, size_type /* j */) {

  1765             return false;

  1766         }

  1767         static

  1768         BOOST_UBLAS_INLINE

  1769         bool other (size_type /* i */, size_type /* j */) {

  1770             return true;

  1771         }

  1772     };

  1774     template <class Z>

  1775     struct basic_lower {

  1776         typedef Z size_type;

  1778         template<class L>

  1779         static

  1780         BOOST_UBLAS_INLINE

  1781         size_type packed_size (L, size_type size1, size_type size2) {

  1782             return L::triangular_size (size1, size2);

  1783         }

  1785         static

  1786         BOOST_UBLAS_INLINE

  1787         bool zero (size_type i, size_type j) {

  1788             return j > i;

  1789         }

  1790         static

  1791         BOOST_UBLAS_INLINE

  1792         bool one (size_type /* i */, size_type /* j */) {

  1793             return false;

  1794         }

  1795         static

  1796         BOOST_UBLAS_INLINE

  1797         bool other (size_type i, size_type j) {

  1798             return j <= i;

  1799         }

  1800         template<class L>

  1801         static

  1802         BOOST_UBLAS_INLINE

  1803         size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {

  1804             return L::lower_element (i, size1, j, size2);

  1805         }

  1807         static

  1808         BOOST_UBLAS_INLINE

  1809         size_type restrict1 (size_type i, size_type j) {

  1810             return (std::max) (i, j);

  1811         }

  1812         static

  1813         BOOST_UBLAS_INLINE

  1814         size_type restrict2 (size_type i, size_type j) {

  1815             return (std::min) (i + 1, j);

  1816         }

  1817         static

  1818         BOOST_UBLAS_INLINE

  1819         size_type mutable_restrict1 (size_type i, size_type j) {

  1820             return (std::max) (i, j);

  1821         }

  1822         static

  1823         BOOST_UBLAS_INLINE

  1824         size_type mutable_restrict2 (size_type i, size_type j) {

  1825             return (std::min) (i + 1, j);

  1826         }

  1827     };

  1828     template <class Z>

  1829     struct basic_upper {

  1830         typedef Z size_type;

  1832         template<class L>

  1833         static

  1834         BOOST_UBLAS_INLINE

  1835         size_type packed_size (L, size_type size1, size_type size2) {

  1836             return L::triangular_size (size1, size2);

  1837         }

  1839         static

  1840         BOOST_UBLAS_INLINE

  1841         bool zero (size_type i, size_type j) {

  1842             return j < i;

  1843         }

  1844         static

  1845         BOOST_UBLAS_INLINE

  1846         bool one (size_type /* i */, size_type /* j */) {

  1847             return false;

  1848         }

  1849         static

  1850         BOOST_UBLAS_INLINE

  1851         bool other (size_type i, size_type j) {

  1852             return j >= i;

  1853         }

  1854         template<class L>

  1855         static

  1856         BOOST_UBLAS_INLINE

  1857         size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {

  1858             return L::upper_element (i, size1, j, size2);

  1859         }

  1861         static

  1862         BOOST_UBLAS_INLINE

  1863         size_type restrict1 (size_type i, size_type j) {

  1864             return (std::min) (i, j + 1);

  1865         }

  1866         static

  1867         BOOST_UBLAS_INLINE

  1868         size_type restrict2 (size_type i, size_type j) {

  1869             return (std::max) (i, j);

  1870         }

  1871         static

  1872         BOOST_UBLAS_INLINE

  1873         size_type mutable_restrict1 (size_type i, size_type j) {

  1874             return (std::min) (i, j + 1);

  1875         }

  1876         static

  1877         BOOST_UBLAS_INLINE

  1878         size_type mutable_restrict2 (size_type i, size_type j) {

  1879             return (std::max) (i, j);

  1880         }

  1881     };

  1882     template <class Z>

  1883     struct basic_unit_lower : public basic_lower<Z> {

  1884         typedef Z size_type;

  1886         template<class L>

  1887         static

  1888         BOOST_UBLAS_INLINE

  1889         size_type packed_size (L, size_type size1, size_type size2) {

  1890             // Zero size strict triangles are bad at this point

  1891             BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());

  1892             return L::triangular_size (size1 - 1, size2 - 1);

  1893         }

  1895         static

  1896         BOOST_UBLAS_INLINE

  1897         bool one (size_type i, size_type j) {

  1898             return j == i;

  1899         }

  1900         static

  1901         BOOST_UBLAS_INLINE

  1902         bool other (size_type i, size_type j) {

  1903             return j < i;

  1904         }

  1905         template<class L>

  1906         static

  1907         BOOST_UBLAS_INLINE

  1908         size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {

  1909             // Zero size strict triangles are bad at this point

  1910             BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());

  1911             return L::lower_element (i, size1 - 1, j, size2 - 1);

  1912         }

  1914         static

  1915         BOOST_UBLAS_INLINE

  1916         size_type mutable_restrict1 (size_type i, size_type j) {

  1917             return (std::max) (i, j);

  1918         }

  1919         static

  1920         BOOST_UBLAS_INLINE

  1921         size_type mutable_restrict2 (size_type i, size_type j) {

  1922             return (std::min) (i, j);

  1923         }

  1924     };

  1925     template <class Z>

  1926     struct basic_unit_upper : public basic_upper<Z> {

  1927         typedef Z size_type;

  1929         template<class L>

  1930         static

  1931         BOOST_UBLAS_INLINE

  1932         size_type packed_size (L, size_type size1, size_type size2) {

  1933             // Zero size strict triangles are bad at this point

  1934             BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());

  1935             return L::triangular_size (size1 - 1, size2 - 1);

  1936         }

  1938         static

  1939         BOOST_UBLAS_INLINE

  1940         bool one (size_type i, size_type j) {

  1941             return j == i;

  1942         }

  1943         static

  1944         BOOST_UBLAS_INLINE

  1945         bool other (size_type i, size_type j) {

  1946             return j > i;

  1947         }

  1948         template<class L>

  1949         static

  1950         BOOST_UBLAS_INLINE

  1951         size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {

  1952             // Zero size strict triangles are bad at this point

  1953             BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());

  1954             return L::upper_element (i, size1 - 1, j, size2 - 1);

  1955         }

  1957         static

  1958         BOOST_UBLAS_INLINE

  1959         size_type mutable_restrict1 (size_type i, size_type j) {

  1960             return (std::min) (i, j);

  1961         }

  1962         static

  1963         BOOST_UBLAS_INLINE

  1964         size_type mutable_restrict2 (size_type i, size_type j) {

  1965             return (std::max) (i, j);

  1966         }

  1967     };

  1968     template <class Z>

  1969     struct basic_strict_lower : public basic_lower<Z> {

  1970         typedef Z size_type;

  1972         template<class L>

  1973         static

  1974         BOOST_UBLAS_INLINE

  1975         size_type packed_size (L, size_type size1, size_type size2) {

  1976             // Zero size strict triangles are bad at this point

  1977             BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());

  1978             return L::triangular_size (size1 - 1, size2 - 1);

  1979         }

  1981         static

  1982         BOOST_UBLAS_INLINE

  1983         bool zero (size_type i, size_type j) {

  1984             return j >= i;

  1985         }

  1986         static

  1987         BOOST_UBLAS_INLINE

  1988         bool one (size_type /*i*/, size_type /*j*/) {

  1989             return false;

  1990         }

  1991         static

  1992         BOOST_UBLAS_INLINE

  1993         bool other (size_type i, size_type j) {

  1994             return j < i;

  1995         }

  1996         template<class L>

  1997         static

  1998         BOOST_UBLAS_INLINE

  1999         size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {

  2000             // Zero size strict triangles are bad at this point

  2001             BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());

  2002             return L::lower_element (i, size1 - 1, j, size2 - 1);

  2003         }

  2005         static

  2006         BOOST_UBLAS_INLINE

  2007         size_type restrict1 (size_type i, size_type j) {

  2008             return (std::max) (i, j);

  2009         }

  2010         static

  2011         BOOST_UBLAS_INLINE

  2012         size_type restrict2 (size_type i, size_type j) {

  2013             return (std::min) (i, j);

  2014         }

  2015         static

  2016         BOOST_UBLAS_INLINE

  2017         size_type mutable_restrict1 (size_type i, size_type j) {

  2018             return (std::max) (i, j);

  2019         }

  2020         static

  2021         BOOST_UBLAS_INLINE

  2022         size_type mutable_restrict2 (size_type i, size_type j) {

  2023             return (std::min) (i, j);

  2024         }

  2025     };

  2026     template <class Z>

  2027     struct basic_strict_upper : public basic_upper<Z> {

  2028         typedef Z size_type;

  2030         template<class L>

  2031         static

  2032         BOOST_UBLAS_INLINE

  2033         size_type packed_size (L, size_type size1, size_type size2) {

  2034             // Zero size strict triangles are bad at this point

  2035             BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());

  2036             return L::triangular_size (size1 - 1, size2 - 1);

  2037         }

  2039         static

  2040         BOOST_UBLAS_INLINE

  2041         bool zero (size_type i, size_type j) {

  2042             return j <= i;

  2043         }

  2044         static

  2045         BOOST_UBLAS_INLINE

  2046         bool one (size_type /*i*/, size_type /*j*/) {

  2047             return false;

  2048         }

  2049         static

  2050         BOOST_UBLAS_INLINE

  2051         bool other (size_type i, size_type j) {

  2052             return j > i;

  2053         }

  2054         template<class L>

  2055         static

  2056         BOOST_UBLAS_INLINE

  2057         size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {

  2058             // Zero size strict triangles are bad at this point

  2059             BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());

  2060             return L::upper_element (i, size1 - 1, j, size2 - 1);

  2061         }

  2063         static

  2064         BOOST_UBLAS_INLINE

  2065         size_type restrict1 (size_type i, size_type j) {

  2066             return (std::min) (i, j);

  2067         }

  2068         static

  2069         BOOST_UBLAS_INLINE

  2070         size_type restrict2 (size_type i, size_type j) {

  2071             return (std::max) (i, j);

  2072         }

  2073         static

  2074         BOOST_UBLAS_INLINE

  2075         size_type mutable_restrict1 (size_type i, size_type j) {

  2076             return (std::min) (i, j);

  2077         }

  2078         static

  2079         BOOST_UBLAS_INLINE

  2080         size_type mutable_restrict2 (size_type i, size_type j) {

  2081             return (std::max) (i, j);

  2082         }

  2083     };

  2085 }}}

  2087 #endif

author	sl
	Tue, 10 Jun 2014 14:32:02 +0200
changeset 1	260cb5ec6c19
permissions	-rw-r--r--