FM / ublas

/usr/include/boost/numeric/ublas

Name	Size	Mode	Actions
detail/	-	0755	rm
experimental/	-	0755	rm
opencl/	-	0755	rm
operation/	-	0755	rm
tensor/	-	0755	rm
traits/	-	0755	rm
assignment.hpp	35756	0644	edit dl rm
banded.hpp	86847	0644	edit dl rm
blas.hpp	18132	0644	edit dl rm
doxydoc.hpp	4087	0644	edit dl rm
exception.hpp	8822	0644	edit dl rm
expression_types.hpp	18467	0644	edit dl rm
functional.hpp	80865	0644	edit dl rm
fwd.hpp	7215	0644	edit dl rm
hermitian.hpp	102000	0644	edit dl rm
io.hpp	13687	0644	edit dl rm
lu.hpp	14044	0644	edit dl rm
matrix.hpp	216096	0644	edit dl rm
matrix_expression.hpp	169092	0644	edit dl rm
matrix_proxy.hpp	201494	0644	edit dl rm
matrix_sparse.hpp	231409	0644	edit dl rm
matrix_vector.hpp	11475	0644	edit dl rm
opencl.hpp	490	0644	edit dl rm
operation.hpp	32110	0644	edit dl rm
operations.hpp	739	0644	edit dl rm
operation_blocked.hpp	13813	0644	edit dl rm
operation_sparse.hpp	8146	0644	edit dl rm
storage.hpp	67562	0644	edit dl rm
storage_sparse.hpp	19762	0644	edit dl rm
symmetric.hpp	88183	0644	edit dl rm
tags.hpp	731	0644	edit dl rm
tensor.hpp	800	0644	edit dl rm
traits.hpp	24952	0644	edit dl rm
triangular.hpp	105304	0644	edit dl rm
vector.hpp	105769	0644	edit dl rm
vector_expression.hpp	65815	0644	edit dl rm
vector_of_vector.hpp	53270	0644	edit dl rm
vector_proxy.hpp	64597	0644	edit dl rm
vector_sparse.hpp	84156	0644	edit dl rm

Edit: /usr/include/boost/numeric/ublas/functional.hpp (80865B)

#ifndef _BOOST_UBLAS_FUNCTIONAL_
#define _BOOST_UBLAS_FUNCTIONAL_

#include <functional>

#include <boost/core/ignore_unused.hpp>

#include <boost/numeric/ublas/traits.hpp>
#ifdef BOOST_UBLAS_USE_DUFF_DEVICE
#include <boost/numeric/ublas/detail/duff.hpp>
#endif
#ifdef BOOST_UBLAS_USE_SIMD
#include <boost/numeric/ublas/detail/raw.hpp>
#else
namespace boost { namespace numeric { namespace ublas { namespace raw {
}}}}
#endif
#ifdef BOOST_UBLAS_HAVE_BINDINGS
#include <boost/numeric/bindings/traits/std_vector.hpp>
#include <boost/numeric/bindings/traits/ublas_vector.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/bindings/atlas/cblas.hpp>
#endif

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

namespace boost { namespace numeric { namespace ublas {

// Scalar functors

// Unary
    template<class T>
    struct scalar_unary_functor {
        typedef T value_type;
        typedef typename type_traits<T>::const_reference argument_type;
        typedef typename type_traits<T>::value_type result_type;
    };

template<class T>
    struct scalar_identity:
        public scalar_unary_functor<T> {
        typedef typename scalar_unary_functor<T>::argument_type argument_type;
        typedef typename scalar_unary_functor<T>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument_type t) {
            return t;
        }
    };
    template<class T>
    struct scalar_negate:
        public scalar_unary_functor<T> {
        typedef typename scalar_unary_functor<T>::argument_type argument_type;
        typedef typename scalar_unary_functor<T>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument_type t) {
            return - t;
        }
    };
    template<class T>
    struct scalar_conj:
        public scalar_unary_functor<T> {
        typedef typename scalar_unary_functor<T>::value_type value_type;
        typedef typename scalar_unary_functor<T>::argument_type argument_type;
        typedef typename scalar_unary_functor<T>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument_type t) {
            return type_traits<value_type>::conj (t);
        }
    };

// Unary returning real
    template<class T>
    struct scalar_real_unary_functor {
        typedef T value_type;
        typedef typename type_traits<T>::const_reference argument_type;
        typedef typename type_traits<T>::real_type result_type;
    };

template<class T>
    struct scalar_real:
        public scalar_real_unary_functor<T> {
        typedef typename scalar_real_unary_functor<T>::value_type value_type;
        typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
        typedef typename scalar_real_unary_functor<T>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument_type t) {
            return type_traits<value_type>::real (t);
        }
    };
    template<class T>
    struct scalar_imag:
        public scalar_real_unary_functor<T> {
        typedef typename scalar_real_unary_functor<T>::value_type value_type;
        typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
        typedef typename scalar_real_unary_functor<T>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument_type t) {
            return type_traits<value_type>::imag (t);
        }
    };

// Binary
    template<class T1, class T2>
    struct scalar_binary_functor {
        typedef typename type_traits<T1>::const_reference argument1_type;
        typedef typename type_traits<T2>::const_reference argument2_type;
        typedef typename promote_traits<T1, T2>::promote_type result_type;
    };

template<class T1, class T2>
    struct scalar_plus:
        public scalar_binary_functor<T1, T2> {
        typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
        typedef typename scalar_binary_functor<T1, T2>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument1_type t1, argument2_type t2) {
            return t1 + t2;
        }
    };
    template<class T1, class T2>
    struct scalar_minus:
        public scalar_binary_functor<T1, T2> {
        typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
        typedef typename scalar_binary_functor<T1, T2>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument1_type t1, argument2_type t2) {
            return t1 - t2;
        }
    };
    template<class T1, class T2>
    struct scalar_multiplies:
        public scalar_binary_functor<T1, T2> {
        typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
        typedef typename scalar_binary_functor<T1, T2>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument1_type t1, argument2_type t2) {
            return t1 * t2;
        }
    };
    template<class T1, class T2>
    struct scalar_divides:
        public scalar_binary_functor<T1, T2> {
        typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
        typedef typename scalar_binary_functor<T1, T2>::result_type result_type;

static BOOST_UBLAS_INLINE
        result_type apply (argument1_type t1, argument2_type t2) {
            return t1 / t2;
        }
    };

template<class T1, class T2>
    struct scalar_binary_assign_functor {
        // ISSUE Remove reference to avoid reference to reference problems
        typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
        typedef typename type_traits<T2>::const_reference argument2_type;
    };

struct assign_tag {};
    struct computed_assign_tag {};

template<class T1, class T2>
    struct scalar_assign:
        public scalar_binary_assign_functor<T1, T2> {
        typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
        static const bool computed ;
#else
        static const bool computed = false ;
#endif

static BOOST_UBLAS_INLINE
        void apply (argument1_type t1, argument2_type t2) {
            t1 = t2;
        }

template<class U1, class U2>
        struct rebind {
            typedef scalar_assign<U1, U2> other;
        };
    };

#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    template<class T1, class T2>
    const bool scalar_assign<T1,T2>::computed = false;
#endif

template<class T1, class T2>
    struct scalar_plus_assign:
        public scalar_binary_assign_functor<T1, T2> {
        typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
        static const bool computed ;
#else
        static const bool computed = true ;
#endif

static BOOST_UBLAS_INLINE
        void apply (argument1_type t1, argument2_type t2) {
            t1 += t2;
        }

template<class U1, class U2>
        struct rebind {
            typedef scalar_plus_assign<U1, U2> other;
        };
    };

#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    template<class T1, class T2>
    const bool scalar_plus_assign<T1,T2>::computed = true;
#endif

template<class T1, class T2>
    struct scalar_minus_assign:
        public scalar_binary_assign_functor<T1, T2> {
        typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
        static const bool computed ;
#else
        static const bool computed = true ;
#endif

static BOOST_UBLAS_INLINE
        void apply (argument1_type t1, argument2_type t2) {
            t1 -= t2;
        }

template<class U1, class U2>
        struct rebind {
            typedef scalar_minus_assign<U1, U2> other;
        };
    };

#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
    template<class T1, class T2>
    const bool scalar_minus_assign<T1,T2>::computed = true;
#endif

template<class T1, class T2>
    struct scalar_multiplies_assign:
        public scalar_binary_assign_functor<T1, T2> {
        typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
        static const bool computed = true;

static BOOST_UBLAS_INLINE
        void apply (argument1_type t1, argument2_type t2) {
            t1 *= t2;
        }

template<class U1, class U2>
        struct rebind {
            typedef scalar_multiplies_assign<U1, U2> other;
        };
    };
    template<class T1, class T2>
    struct scalar_divides_assign:
        public scalar_binary_assign_functor<T1, T2> {
        typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
        static const bool computed ;

static BOOST_UBLAS_INLINE
        void apply (argument1_type t1, argument2_type t2) {
            t1 /= t2;
        }

template<class U1, class U2>
        struct rebind {
            typedef scalar_divides_assign<U1, U2> other;
        };
    };
    template<class T1, class T2>
    const bool scalar_divides_assign<T1,T2>::computed = true;

template<class T1, class T2>
    struct scalar_binary_swap_functor {
        typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
        typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type;
    };

template<class T1, class T2>
    struct scalar_swap:
        public scalar_binary_swap_functor<T1, T2> {
        typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type;
        typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type;

static BOOST_UBLAS_INLINE
        void apply (argument1_type t1, argument2_type t2) {
            std::swap (t1, t2);
        }

template<class U1, class U2>
        struct rebind {
            typedef scalar_swap<U1, U2> other;
        };
    };

// Vector functors

// Unary returning scalar
    template<class V>
    struct vector_scalar_unary_functor {
        typedef typename V::value_type value_type;
        typedef typename V::value_type result_type;
    };

template<class V>
    struct vector_sum: 
        public vector_scalar_unary_functor<V> {
        typedef typename vector_scalar_unary_functor<V>::value_type value_type;
        typedef typename vector_scalar_unary_functor<V>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E> &e) { 
            result_type t = result_type (0);
            typedef typename E::size_type vector_size_type;
            vector_size_type size (e ().size ());
            for (vector_size_type i = 0; i < size; ++ i)
                t += e () (i);
            return t;
        }
        // Dense case
        template<class D, class I>
        static BOOST_UBLAS_INLINE
        result_type apply (D size, I it) { 
            result_type t = result_type (0);
            while (-- size >= 0)
                t += *it, ++ it;
            return t; 
        }
        // Sparse case
        template<class I>
        static BOOST_UBLAS_INLINE
        result_type apply (I it, const I &it_end) {
            result_type t = result_type (0);
            while (it != it_end) 
                t += *it, ++ it;
            return t; 
        }
    };

// Unary returning real scalar 
    template<class V>
    struct vector_scalar_real_unary_functor {
        typedef typename V::value_type value_type;
        typedef typename type_traits<value_type>::real_type real_type;
        typedef real_type result_type;
    };

template<class V>
    struct vector_norm_1:
        public vector_scalar_real_unary_functor<V> {
        typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
        typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
        typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E> &e) {
            real_type t = real_type ();
            typedef typename E::size_type vector_size_type;
            vector_size_type size (e ().size ());
            for (vector_size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::type_abs (e () (i)));
                t += u;
            }
            return t;
        }
        // Dense case
        template<class D, class I>
        static BOOST_UBLAS_INLINE
        result_type apply (D size, I it) {
            real_type t = real_type ();
            while (-- size >= 0) {
                real_type u (type_traits<value_type>::norm_1 (*it));
                t += u;
                ++ it;
            }
            return t;
        }
        // Sparse case
        template<class I>
        static BOOST_UBLAS_INLINE
        result_type apply (I it, const I &it_end) {
            real_type t = real_type ();
            while (it != it_end) {
                real_type u (type_traits<value_type>::norm_1 (*it));
                t += u;
                ++ it;
            }
            return t;
        }
    };
    template<class V>
    struct vector_norm_2:
        public vector_scalar_real_unary_functor<V> {
        typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
        typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
        typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E> &e) {
            typedef typename E::size_type vector_size_type;
            vector_size_type size (e ().size ());
#ifndef BOOST_UBLAS_SCALED_NORM
            real_type t = real_type ();
            for (vector_size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_2 (e () (i)));
                t +=  u * u;
            }
            return static_cast<result_type>(type_traits<real_type>::type_sqrt (t));
#else
            real_type scale = real_type ();
            real_type sum_squares (1);
            for (vector_size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_2 (e () (i)));
                if ( real_type () /* zero */ == u ) continue;
                if (scale < u) {
                    real_type v (scale / u);
                    sum_squares = sum_squares * v * v + real_type (1);
                    scale = u;
                } else {
                    real_type v (u / scale);
                    sum_squares += v * v;
                }
            }
            return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares));
#endif
        }
        // Dense case
        template<class D, class I>
        static BOOST_UBLAS_INLINE
        result_type apply (D size, I it) {
#ifndef BOOST_UBLAS_SCALED_NORM
            real_type t = real_type ();
            while (-- size >= 0) {
                real_type u (type_traits<value_type>::norm_2 (*it));
                t +=  u * u;
                ++ it;
            }
            return static_cast<result_type>(type_traits<real_type>::type_sqrt (t));
#else
            real_type scale = real_type ();
            real_type sum_squares (1);
            while (-- size >= 0) {
                real_type u (type_traits<value_type>::norm_2 (*it));
                if (scale < u) {
                    real_type v (scale / u);
                    sum_squares = sum_squares * v * v + real_type (1);
                    scale = u;
                } else {
                    real_type v (u / scale);
                    sum_squares += v * v;
                }
                ++ it;
            }
            return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares));
#endif
        }
        // Sparse case
        template<class I>
        static BOOST_UBLAS_INLINE
        result_type apply (I it, const I &it_end) {
#ifndef BOOST_UBLAS_SCALED_NORM
            real_type t = real_type ();
            while (it != it_end) {
                real_type u (type_traits<value_type>::norm_2 (*it));
                t +=  u * u;
                ++ it;
            }
            return static_cast<result_type>(type_traits<real_type>::type_sqrt (t));
#else
            real_type scale = real_type ();
            real_type sum_squares (1);
            while (it != it_end) {
                real_type u (type_traits<value_type>::norm_2 (*it));
                if (scale < u) {
                    real_type v (scale / u);
                    sum_squares = sum_squares * v * v + real_type (1);
                    scale = u;
                } else {
                    real_type v (u / scale);
                    sum_squares += v * v;
                }
                ++ it;
            }
            return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares));
#endif
        }
    };

template<class V>
    struct vector_norm_2_square :
        public vector_scalar_real_unary_functor<V> {
        typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
        typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
        typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E> &e) {
            real_type t = real_type ();
            typedef typename E::size_type vector_size_type;
            vector_size_type size (e ().size ());
            for (vector_size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_2 (e () (i)));
                t +=  u * u;
            }
            return t;
        }
        // Dense case
        template<class D, class I>
        static BOOST_UBLAS_INLINE
        result_type apply (D size, I it) {
            real_type t = real_type ();
            while (-- size >= 0) {
                real_type u (type_traits<value_type>::norm_2 (*it));
                t +=  u * u;
                ++ it;
            }
            return t;
        }
        // Sparse case
        template<class I>
        static BOOST_UBLAS_INLINE
        result_type apply (I it, const I &it_end) {
            real_type t = real_type ();
            while (it != it_end) {
                real_type u (type_traits<value_type>::norm_2 (*it));
                t +=  u * u;
                ++ it;
            }
            return t;
        }
    };

template<class V>
    struct vector_norm_inf:
        public vector_scalar_real_unary_functor<V> {
        typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
        typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
        typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E> &e) {
            real_type t = real_type ();
            typedef typename E::size_type vector_size_type;
            vector_size_type size (e ().size ());
            for (vector_size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_inf (e () (i)));
                if (u > t)
                    t = u;
            }
            return t;
        }
        // Dense case
        template<class D, class I>
        static BOOST_UBLAS_INLINE
        result_type apply (D size, I it) {
            real_type t = real_type ();
            while (-- size >= 0) {
                real_type u (type_traits<value_type>::norm_inf (*it));
                if (u > t)
                    t = u;
                ++ it;
            }
            return t;
        }
        // Sparse case
        template<class I>
        static BOOST_UBLAS_INLINE
        result_type apply (I it, const I &it_end) { 
            real_type t = real_type ();
            while (it != it_end) {
                real_type u (type_traits<value_type>::norm_inf (*it));
                if (u > t) 
                    t = u;
                ++ it;
            }
            return t; 
        }
    };

// Unary returning index
    template<class V>
    struct vector_scalar_index_unary_functor {
        typedef typename V::value_type value_type;
        typedef typename type_traits<value_type>::real_type real_type;
        typedef typename V::size_type result_type;
    };

template<class V>
    struct vector_index_norm_inf:
        public vector_scalar_index_unary_functor<V> {
        typedef typename vector_scalar_index_unary_functor<V>::value_type value_type;
        typedef typename vector_scalar_index_unary_functor<V>::real_type real_type;
        typedef typename vector_scalar_index_unary_functor<V>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E> &e) {
            // ISSUE For CBLAS compatibility return 0 index in empty case
            result_type i_norm_inf (0);
            real_type t = real_type ();
            typedef typename E::size_type vector_size_type;
            vector_size_type size (e ().size ());
            for (vector_size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_inf (e () (i)));
                if (u > t) {
                    i_norm_inf = i;
                    t = u;
                }
            }
            return i_norm_inf;
        }
        // Dense case
        template<class D, class I>
        static BOOST_UBLAS_INLINE
        result_type apply (D size, I it) {
            // ISSUE For CBLAS compatibility return 0 index in empty case
            result_type i_norm_inf (0);
            real_type t = real_type ();
            while (-- size >= 0) {
                real_type u (type_traits<value_type>::norm_inf (*it));
                if (u > t) {
                    i_norm_inf = it.index ();
                    t = u;
                }
                ++ it;
            }
            return i_norm_inf;
        }
        // Sparse case
        template<class I>
        static BOOST_UBLAS_INLINE
        result_type apply (I it, const I &it_end) {
            // ISSUE For CBLAS compatibility return 0 index in empty case
            result_type i_norm_inf (0);
            real_type t = real_type ();
            while (it != it_end) {
                real_type u (type_traits<value_type>::norm_inf (*it));
                if (u > t) {
                    i_norm_inf = it.index ();
                    t = u;
                }
                ++ it;
            }
            return i_norm_inf;
        }
    };

// Binary returning scalar
    template<class V1, class V2, class TV>
    struct vector_scalar_binary_functor {
        typedef TV value_type;
        typedef TV result_type;
    };

template<class V1, class V2, class TV>
    struct vector_inner_prod:
        public vector_scalar_binary_functor<V1, V2, TV> {
        typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type;
        typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type;

template<class C1, class C2>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_container<C1> &c1,
                           const vector_container<C2> &c2) {
#ifdef BOOST_UBLAS_USE_SIMD
            using namespace raw;
            typedef typename C1::size_type vector_size_type;
            vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
            const typename V1::value_type *data1 = data_const (c1 ());
            const typename V1::value_type *data2 = data_const (c2 ());
            vector_size_type s1 = stride (c1 ());
            vector_size_type s2 = stride (c2 ());
            result_type t = result_type (0);
            if (s1 == 1 && s2 == 1) {
                for (vector_size_type i = 0; i < size; ++ i)
                    t += data1 [i] * data2 [i];
            } else if (s2 == 1) {
                for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
                    t += data1 [i1] * data2 [i];
            } else if (s1 == 1) {
                for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
                    t += data1 [i] * data2 [i2];
            } else {
                for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
                    t += data1 [i1] * data2 [i2];
            }
            return t;
#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
            return boost::numeric::bindings::atlas::dot (c1 (), c2 ());
#else
            return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2));
#endif
        }
        template<class E1, class E2>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E1> &e1,
                           const vector_expression<E2> &e2) {
            typedef typename E1::size_type vector_size_type;
            vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
            result_type t = result_type (0);
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            for (vector_size_type i = 0; i < size; ++ i)
                t += e1 () (i) * e2 () (i);
#else
            vector_size_type i (0);
            DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));
#endif
            return t;
        }
        // Dense case
        template<class D, class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (D size, I1 it1, I2 it2) {
            result_type t = result_type (0);
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            while (-- size >= 0)
                t += *it1 * *it2, ++ it1, ++ it2;
#else
            DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
#endif
            return t;
        }
        // Packed case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
            result_type t = result_type (0);
            typedef typename I1::difference_type vector_difference_type;
            vector_difference_type it1_size (it1_end - it1);
            vector_difference_type it2_size (it2_end - it2);
            vector_difference_type diff (0);
            if (it1_size > 0 && it2_size > 0)
                diff = it2.index () - it1.index ();
            if (diff != 0) {
                vector_difference_type size = (std::min) (diff, it1_size);
                if (size > 0) {
                    it1 += size;
                    it1_size -= size;
                    diff -= size;
                }
                size = (std::min) (- diff, it2_size);
                if (size > 0) {
                    it2 += size;
                    it2_size -= size;
                    diff += size;
                }
            }
            vector_difference_type size ((std::min) (it1_size, it2_size));
            while (-- size >= 0)
                t += *it1 * *it2, ++ it1, ++ it2;
            return t;
        }
        // Sparse case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
            result_type t = result_type (0);
            if (it1 != it1_end && it2 != it2_end) {
                for (;;) {
                    if (it1.index () == it2.index ()) {
                        t += *it1 * *it2, ++ it1, ++ it2;
                        if (it1 == it1_end || it2 == it2_end)
                            break;
                    } else if (it1.index () < it2.index ()) {
                        increment (it1, it1_end, it2.index () - it1.index ());
                        if (it1 == it1_end)
                            break;
                    } else if (it1.index () > it2.index ()) {
                        increment (it2, it2_end, it1.index () - it2.index ());
                        if (it2 == it2_end)
                            break;
                    }
                }
            }
            return t;
        }
    };

// Matrix functors

// Binary returning vector
    template<class M1, class M2, class TV>
    struct matrix_vector_binary_functor {
        typedef typename M1::size_type size_type;
        typedef typename M1::difference_type difference_type;
        typedef TV value_type;
        typedef TV result_type;
    };

template<class M1, class M2, class TV>
    struct matrix_vector_prod1:
        public matrix_vector_binary_functor<M1, M2, TV> {
        typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
        typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
        typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
        typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;

template<class C1, class C2>
        static BOOST_UBLAS_INLINE
        result_type apply (const matrix_container<C1> &c1,
                           const vector_container<C2> &c2,
                           size_type i) {
#ifdef BOOST_UBLAS_USE_SIMD
            using namespace raw;
            size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ());
            const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
            const typename M2::value_type *data2 = data_const (c2 ());
            size_type s1 = stride2 (c1 ());
            size_type s2 = stride (c2 ());
            result_type t = result_type (0);
            if (s1 == 1 && s2 == 1) {
                for (size_type j = 0; j < size; ++ j)
                    t += data1 [j] * data2 [j];
            } else if (s2 == 1) {
                for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
                    t += data1 [j1] * data2 [j];
            } else if (s1 == 1) {
                for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
                    t += data1 [j] * data2 [j2];
            } else {
                for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
                    t += data1 [j1] * data2 [j2];
            }
            return t;
#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
            return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ());
#else
            return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i));
#endif
        }
        template<class E1, class E2>
        static BOOST_UBLAS_INLINE
        result_type apply (const matrix_expression<E1> &e1,
                           const vector_expression<E2> &e2,
                           size_type i) {
            size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
            result_type t = result_type (0);
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            for (size_type j = 0; j < size; ++ j)
                t += e1 () (i, j) * e2 () (j);
#else
            size_type j (0);
            DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j));
#endif
            return t;
        }
        // Dense case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (difference_type size, I1 it1, I2 it2) {
            result_type t = result_type (0);
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            while (-- size >= 0)
                t += *it1 * *it2, ++ it1, ++ it2;
#else
            DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
#endif
            return t;
        }
        // Packed case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
            result_type t = result_type (0);
            difference_type it1_size (it1_end - it1);
            difference_type it2_size (it2_end - it2);
            difference_type diff (0);
            if (it1_size > 0 && it2_size > 0)
                diff = it2.index () - it1.index2 ();
            if (diff != 0) {
                difference_type size = (std::min) (diff, it1_size);
                if (size > 0) {
                    it1 += size;
                    it1_size -= size;
                    diff -= size;
                }
                size = (std::min) (- diff, it2_size);
                if (size > 0) {
                    it2 += size;
                    it2_size -= size;
                    diff += size;
                }
            }
            difference_type size ((std::min) (it1_size, it2_size));
            while (-- size >= 0)
                t += *it1 * *it2, ++ it1, ++ it2;
            return t;
        }
        // Sparse case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
                           sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
            result_type t = result_type (0);
            if (it1 != it1_end && it2 != it2_end) {
                size_type it1_index = it1.index2 (), it2_index = it2.index ();
                for (;;) {
                    difference_type compare = it1_index - it2_index;
                    if (compare == 0) {
                        t += *it1 * *it2, ++ it1, ++ it2;
                        if (it1 != it1_end && it2 != it2_end) {
                            it1_index = it1.index2 ();
                            it2_index = it2.index ();
                        } else
                            break;
                    } else if (compare < 0) {
                        increment (it1, it1_end, - compare);
                        if (it1 != it1_end)
                            it1_index = it1.index2 ();
                        else
                            break;
                    } else if (compare > 0) {
                        increment (it2, it2_end, compare);
                        if (it2 != it2_end)
                            it2_index = it2.index ();
                        else
                            break;
                    }
                }
            }
            return t;
        }
        // Sparse packed case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
                           sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
            result_type t = result_type (0);
            while (it1 != it1_end) {
                t += *it1 * it2 () (it1.index2 ());
                ++ it1;
            }
            return t;
        }
        // Packed sparse case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
                           packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
            result_type t = result_type (0);
            while (it2 != it2_end) {
                t += it1 () (it1.index1 (), it2.index ()) * *it2;
                ++ it2;
            }
            return t;
        }
        // Another dispatcher
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
                           sparse_bidirectional_iterator_tag) {
            typedef typename I1::iterator_category iterator1_category;
            typedef typename I2::iterator_category iterator2_category;
            return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
        }
    };

template<class M1, class M2, class TV>
    struct matrix_vector_prod2:
        public matrix_vector_binary_functor<M1, M2, TV> {
        typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
        typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
        typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
        typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;

template<class C1, class C2>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_container<C1> &c1,
                           const matrix_container<C2> &c2,
                           size_type i) {
#ifdef BOOST_UBLAS_USE_SIMD
            using namespace raw;
            size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ());
            const typename M1::value_type *data1 = data_const (c1 ());
            const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ());
            size_type s1 = stride (c1 ());
            size_type s2 = stride1 (c2 ());
            result_type t = result_type (0);
            if (s1 == 1 && s2 == 1) {
                for (size_type j = 0; j < size; ++ j)
                    t += data1 [j] * data2 [j];
            } else if (s2 == 1) {
                for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
                    t += data1 [j1] * data2 [j];
            } else if (s1 == 1) {
                for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
                    t += data1 [j] * data2 [j2];
            } else {
                for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
                    t += data1 [j1] * data2 [j2];
            }
            return t;
#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
            return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i));
#else
            return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
#endif
        }
        template<class E1, class E2>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E1> &e1,
                           const matrix_expression<E2> &e2,
                           size_type i) {
            size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
            result_type t = result_type (0);
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            for (size_type j = 0; j < size; ++ j)
                t += e1 () (j) * e2 () (j, i);
#else
            size_type j (0);
            DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j));
#endif
            return t;
        }
        // Dense case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (difference_type size, I1 it1, I2 it2) {
            result_type t = result_type (0);
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            while (-- size >= 0)
                t += *it1 * *it2, ++ it1, ++ it2;
#else
            DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
#endif
            return t;
        }
        // Packed case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
            result_type t = result_type (0);
            difference_type it1_size (it1_end - it1);
            difference_type it2_size (it2_end - it2);
            difference_type diff (0);
            if (it1_size > 0 && it2_size > 0)
                diff = it2.index1 () - it1.index ();
            if (diff != 0) {
                difference_type size = (std::min) (diff, it1_size);
                if (size > 0) {
                    it1 += size;
                    it1_size -= size;
                    diff -= size;
                }
                size = (std::min) (- diff, it2_size);
                if (size > 0) {
                    it2 += size;
                    it2_size -= size;
                    diff += size;
                }
            }
            difference_type size ((std::min) (it1_size, it2_size));
            while (-- size >= 0)
                t += *it1 * *it2, ++ it1, ++ it2;
            return t;
        }
        // Sparse case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
                           sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
            result_type t = result_type (0);
            if (it1 != it1_end && it2 != it2_end) {
                size_type it1_index = it1.index (), it2_index = it2.index1 ();
                for (;;) {
                    difference_type compare = it1_index - it2_index;
                    if (compare == 0) {
                        t += *it1 * *it2, ++ it1, ++ it2;
                        if (it1 != it1_end && it2 != it2_end) {
                            it1_index = it1.index ();
                            it2_index = it2.index1 ();
                        } else
                            break;
                    } else if (compare < 0) {
                        increment (it1, it1_end, - compare);
                        if (it1 != it1_end)
                            it1_index = it1.index ();
                        else
                            break;
                    } else if (compare > 0) {
                        increment (it2, it2_end, compare);
                        if (it2 != it2_end)
                            it2_index = it2.index1 ();
                        else
                            break;
                    }
                }
            }
            return t;
        }
        // Packed sparse case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
                           packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
            result_type t = result_type (0);
            while (it2 != it2_end) {
                t += it1 () (it2.index1 ()) * *it2;
                ++ it2;
            }
            return t;
        }
        // Sparse packed case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
                           sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
            result_type t = result_type (0);
            while (it1 != it1_end) {
                t += *it1 * it2 () (it1.index (), it2.index2 ());
                ++ it1;
            }
            return t;
        }
        // Another dispatcher
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
                           sparse_bidirectional_iterator_tag) {
            typedef typename I1::iterator_category iterator1_category;
            typedef typename I2::iterator_category iterator2_category;
            return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
        }
    };

// Binary returning matrix
    template<class M1, class M2, class TV>
    struct matrix_matrix_binary_functor {
        typedef typename M1::size_type size_type;
        typedef typename M1::difference_type difference_type;
        typedef TV value_type;
        typedef TV result_type;
    };

template<class M1, class M2, class TV>
    struct matrix_matrix_prod:
        public matrix_matrix_binary_functor<M1, M2, TV> {
        typedef typename matrix_matrix_binary_functor<M1, M2, TV>::size_type size_type;
        typedef typename matrix_matrix_binary_functor<M1, M2, TV>::difference_type difference_type;
        typedef typename matrix_matrix_binary_functor<M1, M2, TV>::value_type value_type;
        typedef typename matrix_matrix_binary_functor<M1, M2, TV>::result_type result_type;

template<class C1, class C2>
        static BOOST_UBLAS_INLINE
        result_type apply (const matrix_container<C1> &c1,
                           const matrix_container<C2> &c2,
                           size_type i, size_type j) {
#ifdef BOOST_UBLAS_USE_SIMD
            using namespace raw;
            size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ());
            const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
            const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ());
            size_type s1 = stride2 (c1 ());
            size_type s2 = stride1 (c2 ());
            result_type t = result_type (0);
            if (s1 == 1 && s2 == 1) {
                for (size_type k = 0; k < size; ++ k)
                    t += data1 [k] * data2 [k];
            } else if (s2 == 1) {
                for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1)
                    t += data1 [k1] * data2 [k];
            } else if (s1 == 1) {
                for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2)
                    t += data1 [k] * data2 [k2];
            } else {
                for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2)
                    t += data1 [k1] * data2 [k2];
            }
            return t;
#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
            return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j));
#else
            boost::ignore_unused(j);
            return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
#endif
        }
        template<class E1, class E2>
        static BOOST_UBLAS_INLINE
        result_type apply (const matrix_expression<E1> &e1,
                           const matrix_expression<E2> &e2,
                           size_type i, size_type j) {
            size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
            result_type t = result_type (0);
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            for (size_type k = 0; k < size; ++ k)
                t += e1 () (i, k) * e2 () (k, j);
#else
            size_type k (0);
            DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k));
#endif
            return t;
        }
        // Dense case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (difference_type size, I1 it1, I2 it2) {
            result_type t = result_type (0);
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            while (-- size >= 0)
                t += *it1 * *it2, ++ it1, ++ it2;
#else
            DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
#endif
            return t;
        }
        // Packed case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) {
            result_type t = result_type (0);
            difference_type it1_size (it1_end - it1);
            difference_type it2_size (it2_end - it2);
            difference_type diff (0);
            if (it1_size > 0 && it2_size > 0)
                diff = it2.index1 () - it1.index2 ();
            if (diff != 0) {
                difference_type size = (std::min) (diff, it1_size);
                if (size > 0) {
                    it1 += size;
                    it1_size -= size;
                    diff -= size;
                }
                size = (std::min) (- diff, it2_size);
                if (size > 0) {
                    it2 += size;
                    it2_size -= size;
                    diff += size;
                }
            }
            difference_type size ((std::min) (it1_size, it2_size));
            while (-- size >= 0)
                t += *it1 * *it2, ++ it1, ++ it2;
            return t;
        }
        // Sparse case
        template<class I1, class I2>
        static BOOST_UBLAS_INLINE
        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
            result_type t = result_type (0);
            if (it1 != it1_end && it2 != it2_end) {
                size_type it1_index = it1.index2 (), it2_index = it2.index1 ();
                for (;;) {
                    difference_type compare = difference_type (it1_index - it2_index);
                    if (compare == 0) {
                        t += *it1 * *it2, ++ it1, ++ it2;
                        if (it1 != it1_end && it2 != it2_end) {
                            it1_index = it1.index2 ();
                            it2_index = it2.index1 ();
                        } else
                            break;
                    } else if (compare < 0) {
                        increment (it1, it1_end, - compare);
                        if (it1 != it1_end)
                            it1_index = it1.index2 ();
                        else
                            break;
                    } else if (compare > 0) {
                        increment (it2, it2_end, compare);
                        if (it2 != it2_end)
                            it2_index = it2.index1 ();
                        else
                            break;
                    }
                }
            }
            return t;
        }
    };

// Unary returning scalar norm
    template<class M>
    struct matrix_scalar_real_unary_functor {
        typedef typename M::value_type value_type;
        typedef typename type_traits<value_type>::real_type real_type;
        typedef real_type result_type;
    };

template<class M>
    struct matrix_norm_1:
        public matrix_scalar_real_unary_functor<M> {
        typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
        typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
        typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const matrix_expression<E> &e) {
            real_type t = real_type ();
            typedef typename E::size_type matrix_size_type;
            matrix_size_type size2 (e ().size2 ());
            for (matrix_size_type j = 0; j < size2; ++ j) {
                real_type u = real_type ();
                matrix_size_type size1 (e ().size1 ());
                for (matrix_size_type i = 0; i < size1; ++ i) {
                    real_type v (type_traits<value_type>::norm_1 (e () (i, j)));
                    u += v;
                }
                if (u > t)
                    t = u;
            }
            return t; 
        }
    };

template<class M>
    struct matrix_norm_frobenius:
        public matrix_scalar_real_unary_functor<M> {
        typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
        typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
        typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const matrix_expression<E> &e) { 
            real_type t = real_type ();
            typedef typename E::size_type matrix_size_type;
            matrix_size_type size1 (e ().size1 ());
            for (matrix_size_type i = 0; i < size1; ++ i) {
                matrix_size_type size2 (e ().size2 ());
                for (matrix_size_type j = 0; j < size2; ++ j) {
                    real_type u (type_traits<value_type>::norm_2 (e () (i, j)));
                    t +=  u * u;
                }
            }
            return type_traits<real_type>::type_sqrt (t); 
        }
    };

template<class M>
    struct matrix_norm_inf: 
        public matrix_scalar_real_unary_functor<M> {
        typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
        typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
        typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;

template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const matrix_expression<E> &e) {
            real_type t = real_type ();
            typedef typename E::size_type matrix_size_type;
            matrix_size_type size1 (e ().size1 ());
            for (matrix_size_type i = 0; i < size1; ++ i) {
                real_type u = real_type ();
                matrix_size_type size2 (e ().size2 ());
                for (matrix_size_type j = 0; j < size2; ++ j) {
                    real_type v (type_traits<value_type>::norm_inf (e () (i, j)));
                    u += v;
                }
                if (u > t) 
                    t = u;  
            }
            return t; 
        }
    };

// forward declaration
    template <class Z, class D> struct basic_column_major;

// This functor defines storage layout and it's properties
    // matrix (i,j) -> storage [i * size_i + j]
    template <class Z, class D>
    struct basic_row_major {
        typedef Z size_type;
        typedef D difference_type;
        typedef row_major_tag orientation_category;
        typedef basic_column_major<Z,D> transposed_layout;

static
        BOOST_UBLAS_INLINE
        size_type storage_size (size_type size_i, size_type size_j) {
            // Guard against size_type overflow
            BOOST_UBLAS_CHECK (size_j == 0 || size_i <= (std::numeric_limits<size_type>::max) () / size_j, bad_size ());
            return size_i * size_j;
        }

// Indexing conversion to storage element
        static
        BOOST_UBLAS_INLINE
        size_type element (size_type i, size_type size_i, size_type j, size_type size_j) {
            BOOST_UBLAS_CHECK (i < size_i, bad_index ());
            BOOST_UBLAS_CHECK (j < size_j, bad_index ());
            detail::ignore_unused_variable_warning(size_i);
            // Guard against size_type overflow
            BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ());
            return i * size_j + j;
        }
        static
        BOOST_UBLAS_INLINE
        size_type address (size_type i, size_type size_i, size_type j, size_type size_j) {
            BOOST_UBLAS_CHECK (i <= size_i, bad_index ());
            BOOST_UBLAS_CHECK (j <= size_j, bad_index ());
            // Guard against size_type overflow - address may be size_j past end of storage
            BOOST_UBLAS_CHECK (size_j == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ());
            detail::ignore_unused_variable_warning(size_i);
            return i * size_j + j;
        }

// Storage element to index conversion
        static
        BOOST_UBLAS_INLINE
        difference_type distance_i (difference_type k, size_type /* size_i */, size_type size_j) {
            return size_j != 0 ? k / size_j : 0;
        }
        static
        BOOST_UBLAS_INLINE
        difference_type distance_j (difference_type k, size_type /* size_i */, size_type /* size_j */) {
            return k;
        }
        static
        BOOST_UBLAS_INLINE
        size_type index_i (difference_type k, size_type /* size_i */, size_type size_j) {
            return size_j != 0 ? k / size_j : 0;
        }
        static
        BOOST_UBLAS_INLINE
        size_type index_j (difference_type k, size_type /* size_i */, size_type size_j) {
            return size_j != 0 ? k % size_j : 0;
        }
        static
        BOOST_UBLAS_INLINE
        bool fast_i () {
            return false;
        }
        static
        BOOST_UBLAS_INLINE
        bool fast_j () {
            return true;
        }

// Iterating storage elements
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void increment_i (I &it, size_type /* size_i */, size_type size_j) {
            it += size_j;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void increment_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) {
            it += n * size_j;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void decrement_i (I &it, size_type /* size_i */, size_type size_j) {
            it -= size_j;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) {
            it -= n * size_j;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void increment_j (I &it, size_type /* size_i */, size_type /* size_j */) {
            ++ it;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void increment_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
            it += n;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void decrement_j (I &it, size_type /* size_i */, size_type /* size_j */) {
            -- it;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void decrement_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
            it -= n;
        }

// Triangular access
        static
        BOOST_UBLAS_INLINE
        size_type triangular_size (size_type size_i, size_type size_j) {
            size_type size = (std::max) (size_i, size_j);
            // Guard against size_type overflow - simplified
            BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
            return ((size + 1) * size) / 2;
        }
        static
        BOOST_UBLAS_INLINE
        size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) {
            BOOST_UBLAS_CHECK (i < size_i, bad_index ());
            BOOST_UBLAS_CHECK (j < size_j, bad_index ());
            BOOST_UBLAS_CHECK (i >= j, bad_index ());
            detail::ignore_unused_variable_warning(size_i);
            detail::ignore_unused_variable_warning(size_j);
            // FIXME size_type overflow
            // sigma_i (i + 1) = (i + 1) * i / 2
            // i = 0 1 2 3, sigma = 0 1 3 6
            return ((i + 1) * i) / 2 + j;
        }
        static
        BOOST_UBLAS_INLINE
        size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) {
            BOOST_UBLAS_CHECK (i < size_i, bad_index ());
            BOOST_UBLAS_CHECK (j < size_j, bad_index ());
            BOOST_UBLAS_CHECK (i <= j, bad_index ());
            // FIXME size_type overflow
            // sigma_i (size - i) = size * i - i * (i - 1) / 2
            // i = 0 1 2 3, sigma = 0 4 7 9
            return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i;
        }

// Major and minor indices
        static
        BOOST_UBLAS_INLINE
        size_type index_M (size_type index1, size_type /* index2 */) {
            return index1;
        }
        static
        BOOST_UBLAS_INLINE
        size_type index_m (size_type /* index1 */, size_type index2) {
            return index2;
        }
        static
        BOOST_UBLAS_INLINE
        size_type size_M (size_type size_i, size_type /* size_j */) {
            return size_i;
        }
        static
        BOOST_UBLAS_INLINE
        size_type size_m (size_type /* size_i */, size_type size_j) {
            return size_j;
        }
    };

// This functor defines storage layout and it's properties
    // matrix (i,j) -> storage [i + j * size_i]
    template <class Z, class D>
    struct basic_column_major {
        typedef Z size_type;
        typedef D difference_type;
        typedef column_major_tag orientation_category;
        typedef basic_row_major<Z,D> transposed_layout;

static
        BOOST_UBLAS_INLINE
        size_type storage_size (size_type size_i, size_type size_j) {
            // Guard against size_type overflow
            BOOST_UBLAS_CHECK (size_i == 0 || size_j <= (std::numeric_limits<size_type>::max) () / size_i, bad_size ());
            return size_i * size_j;
        }

// Indexing conversion to storage element
        static
        BOOST_UBLAS_INLINE
        size_type element (size_type i, size_type size_i, size_type j, size_type size_j) {
            BOOST_UBLAS_CHECK (i < size_i, bad_index ());
            BOOST_UBLAS_CHECK (j < size_j, bad_index ());
            detail::ignore_unused_variable_warning(size_j);
            // Guard against size_type overflow
            BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ());
            return i + j * size_i;
        }
        static
        BOOST_UBLAS_INLINE
        size_type address (size_type i, size_type size_i, size_type j, size_type size_j) {
            BOOST_UBLAS_CHECK (i <= size_i, bad_index ());
            BOOST_UBLAS_CHECK (j <= size_j, bad_index ());
            detail::ignore_unused_variable_warning(size_j);
            // Guard against size_type overflow - address may be size_i past end of storage
            BOOST_UBLAS_CHECK (size_i == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ());
            return i + j * size_i;
        }

// Storage element to index conversion
        static
        BOOST_UBLAS_INLINE
        difference_type distance_i (difference_type k, size_type /* size_i */, size_type /* size_j */) {
            return k;
        }
        static
        BOOST_UBLAS_INLINE
        difference_type distance_j (difference_type k, size_type size_i, size_type /* size_j */) {
            return size_i != 0 ? k / size_i : 0;
        }
        static
        BOOST_UBLAS_INLINE
        size_type index_i (difference_type k, size_type size_i, size_type /* size_j */) {
            return size_i != 0 ? k % size_i : 0;
        }
        static
        BOOST_UBLAS_INLINE
        size_type index_j (difference_type k, size_type size_i, size_type /* size_j */) {
            return size_i != 0 ? k / size_i : 0;
        }
        static
        BOOST_UBLAS_INLINE
        bool fast_i () {
            return true;
        }
        static
        BOOST_UBLAS_INLINE
        bool fast_j () {
            return false;
        }

// Iterating
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void increment_i (I &it, size_type /* size_i */, size_type /* size_j */) {
            ++ it;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void increment_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
            it += n;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void decrement_i (I &it, size_type /* size_i */, size_type /* size_j */) {
            -- it;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
            it -= n;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void increment_j (I &it, size_type size_i, size_type /* size_j */) {
            it += size_i;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void increment_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) {
            it += n * size_i;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void decrement_j (I &it, size_type size_i, size_type /* size_j */) {
            it -= size_i;
        }
        template<class I>
        static
        BOOST_UBLAS_INLINE
        void decrement_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) {
            it -= n* size_i;
        }

// Triangular access
        static
        BOOST_UBLAS_INLINE
        size_type triangular_size (size_type size_i, size_type size_j) {
            size_type size = (std::max) (size_i, size_j);
            // Guard against size_type overflow - simplified
            BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
            return ((size + 1) * size) / 2;
        }
        static
        BOOST_UBLAS_INLINE
        size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) {
            BOOST_UBLAS_CHECK (i < size_i, bad_index ());
            BOOST_UBLAS_CHECK (j < size_j, bad_index ());
            BOOST_UBLAS_CHECK (i >= j, bad_index ());
            // FIXME size_type overflow
            // sigma_j (size - j) = size * j - j * (j - 1) / 2
            // j = 0 1 2 3, sigma = 0 4 7 9
            return i - j + (j * (2 * (std::max) (size_i, size_j) - j + 1)) / 2;
        }
        static
        BOOST_UBLAS_INLINE
        size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) {
            BOOST_UBLAS_CHECK (i < size_i, bad_index ());
            BOOST_UBLAS_CHECK (j < size_j, bad_index ());
            BOOST_UBLAS_CHECK (i <= j, bad_index ());
            boost::ignore_unused(size_i, size_j);
            // FIXME size_type overflow
            // sigma_j (j + 1) = (j + 1) * j / 2
            // j = 0 1 2 3, sigma = 0 1 3 6
            return i + ((j + 1) * j) / 2;
        }

// Major and minor indices
        static
        BOOST_UBLAS_INLINE
        size_type index_M (size_type /* index1 */, size_type index2) {
            return index2;
        }
        static
        BOOST_UBLAS_INLINE
        size_type index_m (size_type index1, size_type /* index2 */) {
            return index1;
        }
        static
        BOOST_UBLAS_INLINE
        size_type size_M (size_type /* size_i */, size_type size_j) {
            return size_j;
        }
        static
        BOOST_UBLAS_INLINE
        size_type size_m (size_type size_i, size_type /* size_j */) {
            return size_i;
        }
    };

template <class Z>
    struct basic_full {
        typedef Z size_type;

template<class L>
        static
        BOOST_UBLAS_INLINE
        size_type packed_size (L, size_type size_i, size_type size_j) {
            return L::storage_size (size_i, size_j);
        }

static
        BOOST_UBLAS_INLINE
        bool zero (size_type /* i */, size_type /* j */) {
            return false;
        }
        static
        BOOST_UBLAS_INLINE
        bool one (size_type /* i */, size_type /* j */) {
            return false;
        }
        static
        BOOST_UBLAS_INLINE
        bool other (size_type /* i */, size_type /* j */) {
            return true;
        }
        // FIXME: this should not be used at all
        static
        BOOST_UBLAS_INLINE
        size_type restrict1 (size_type i, size_type /* j */) {
            return i;
        }
        static
        BOOST_UBLAS_INLINE
        size_type restrict2 (size_type /* i */, size_type j) {
            return j;
        }
        static
        BOOST_UBLAS_INLINE
        size_type mutable_restrict1 (size_type i, size_type /* j */) {
            return i;
        }
        static
        BOOST_UBLAS_INLINE
        size_type mutable_restrict2 (size_type /* i */, size_type j) {
            return j;
        }
    };

namespace detail {
        template < class L >
        struct transposed_structure {
            typedef typename L::size_type size_type;

template<class LAYOUT>
            static
            BOOST_UBLAS_INLINE
            size_type packed_size (LAYOUT l, size_type size_i, size_type size_j) {
                return L::packed_size(l, size_j, size_i);
            }

static
            BOOST_UBLAS_INLINE
            bool zero (size_type i, size_type j) {
                return L::zero(j, i);
            }
            static
            BOOST_UBLAS_INLINE
            bool one (size_type i, size_type j) {
                return L::one(j, i);
            }
            static
            BOOST_UBLAS_INLINE
            bool other (size_type i, size_type j) {
                return L::other(j, i);
            }
            template<class LAYOUT>
            static
            BOOST_UBLAS_INLINE
            size_type element (LAYOUT /* l */, size_type i, size_type size_i, size_type j, size_type size_j) {
                return L::element(typename LAYOUT::transposed_layout(), j, size_j, i, size_i);
            }

static
            BOOST_UBLAS_INLINE
            size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
                return L::restrict2(j, i, size2, size1);
            }
            static
            BOOST_UBLAS_INLINE
            size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
                return L::restrict1(j, i, size2, size1);
            }
            static
            BOOST_UBLAS_INLINE
            size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
                return L::mutable_restrict2(j, i, size2, size1);
            }
            static
            BOOST_UBLAS_INLINE
            size_type mutable_restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
                return L::mutable_restrict1(j, i, size2, size1);
            }

static
            BOOST_UBLAS_INLINE
            size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
                return L::global_restrict2(index2, size2, index1, size1);
            }
            static
            BOOST_UBLAS_INLINE
            size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
                return L::global_restrict1(index2, size2, index1, size1);
            }
            static
            BOOST_UBLAS_INLINE
            size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
                return L::global_mutable_restrict2(index2, size2, index1, size1);
            }
            static
            BOOST_UBLAS_INLINE
            size_type global_mutable_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
                return L::global_mutable_restrict1(index2, size2, index1, size1);
            }
        };
    }

template <class Z>
    struct basic_lower {
        typedef Z size_type;
        typedef lower_tag triangular_type;

template<class L>
        static
        BOOST_UBLAS_INLINE
        size_type packed_size (L, size_type size_i, size_type size_j) {
            return L::triangular_size (size_i, size_j);
        }

static
        BOOST_UBLAS_INLINE
        bool zero (size_type i, size_type j) {
            return j > i;
        }
        static
        BOOST_UBLAS_INLINE
        bool one (size_type /* i */, size_type /* j */) {
            return false;
        }
        static
        BOOST_UBLAS_INLINE
        bool other (size_type i, size_type j) {
            return j <= i;
        }
        template<class L>
        static
        BOOST_UBLAS_INLINE
        size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
            return L::lower_element (i, size_i, j, size_j);
        }

// return nearest valid index in column j
        static
        BOOST_UBLAS_INLINE
        size_type restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
            return (std::max)(j, (std::min) (size1, i));
        }
        // return nearest valid index in row i
        static
        BOOST_UBLAS_INLINE
        size_type restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
            return (std::max)(size_type(0), (std::min) (i+1, j));
        }
        // return nearest valid mutable index in column j
        static
        BOOST_UBLAS_INLINE
        size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
            return (std::max)(j, (std::min) (size1, i));
        }
        // return nearest valid mutable index in row i
        static
        BOOST_UBLAS_INLINE
        size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
            return (std::max)(size_type(0), (std::min) (i+1, j));
        }

// return an index between the first and (1+last) filled row
        static
        BOOST_UBLAS_INLINE
        size_type global_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
            return (std::max)(size_type(0), (std::min)(size1, index1) );
        }
        // return an index between the first and (1+last) filled column
        static
        BOOST_UBLAS_INLINE
        size_type global_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
            return (std::max)(size_type(0), (std::min)(size2, index2) );
        }

// return an index between the first and (1+last) filled mutable row
        static
        BOOST_UBLAS_INLINE
        size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
            return (std::max)(size_type(0), (std::min)(size1, index1) );
        }
        // return an index between the first and (1+last) filled mutable column
        static
        BOOST_UBLAS_INLINE
        size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
            return (std::max)(size_type(0), (std::min)(size2, index2) );
        }
    };

// the first row only contains a single 1. Thus it is not stored.
    template <class Z>
    struct basic_unit_lower : public basic_lower<Z> {
        typedef Z size_type;
        typedef unit_lower_tag triangular_type;

template<class L>
        static
        BOOST_UBLAS_INLINE
        size_type packed_size (L, size_type size_i, size_type size_j) {
            // Zero size strict triangles are bad at this point
            BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ());
            return L::triangular_size (size_i - 1, size_j - 1);
        }

static
        BOOST_UBLAS_INLINE
        bool one (size_type i, size_type j) {
            return j == i;
        }
        static
        BOOST_UBLAS_INLINE
        bool other (size_type i, size_type j) {
            return j < i;
        }
        template<class L>
        static
        BOOST_UBLAS_INLINE
        size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
            // Zero size strict triangles are bad at this point
            BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ());
            return L::lower_element (i-1, size_i - 1, j, size_j - 1);
        }

static
        BOOST_UBLAS_INLINE
        size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
            return (std::max)(j+1, (std::min) (size1, i));
        }
        static
        BOOST_UBLAS_INLINE
        size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
            return (std::max)(size_type(0), (std::min) (i, j));
        }

// return an index between the first and (1+last) filled mutable row
        static
        BOOST_UBLAS_INLINE
        size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
            return (std::max)(size_type(1), (std::min)(size1, index1) );
        }
        // return an index between the first and (1+last) filled mutable column
        static
        BOOST_UBLAS_INLINE
        size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
            BOOST_UBLAS_CHECK( size2 >= 1 , external_logic() );
            return (std::max)(size_type(0), (std::min)(size2-1, index2) );
        }
    };

// the first row only contains no element. Thus it is not stored.
    template <class Z>
    struct basic_strict_lower : public basic_unit_lower<Z> {
        typedef Z size_type;
        typedef strict_lower_tag triangular_type;

static
        BOOST_UBLAS_INLINE
        bool zero (size_type i, size_type j) {
            return j >= i;
        }
        static
        BOOST_UBLAS_INLINE
        bool one (size_type /*i*/, size_type /*j*/) {
            return false;
        }
        static
        BOOST_UBLAS_INLINE
        bool other (size_type i, size_type j) {
            return j < i;
        }
        template<class L>
        static
        BOOST_UBLAS_INLINE
        size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
            // Zero size strict triangles are bad at this point
            BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ());
            return L::lower_element (i-1, size_i - 1, j, size_j - 1);
        }

static
        BOOST_UBLAS_INLINE
        size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
            return basic_unit_lower<Z>::mutable_restrict1(i, j, size1, size2);
        }
        static
        BOOST_UBLAS_INLINE
        size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
            return basic_unit_lower<Z>::mutable_restrict2(i, j, size1, size2);
        }

// return an index between the first and (1+last) filled row
        static
        BOOST_UBLAS_INLINE
        size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
            return basic_unit_lower<Z>::global_mutable_restrict1(index1, size1, index2, size2);
        }
        // return an index between the first and (1+last) filled column
        static
        BOOST_UBLAS_INLINE
        size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
            return basic_unit_lower<Z>::global_mutable_restrict2(index1, size1, index2, size2);
        }
    };

template <class Z>
    struct basic_upper : public detail::transposed_structure<basic_lower<Z> >
    { 
        typedef upper_tag triangular_type;
    };

template <class Z>
    struct basic_unit_upper : public detail::transposed_structure<basic_unit_lower<Z> >
    { 
        typedef unit_upper_tag triangular_type;
    };

template <class Z>
    struct basic_strict_upper : public detail::transposed_structure<basic_strict_lower<Z> >
    { 
        typedef strict_upper_tag triangular_type;
    };

}}}

#endif

// // Copyright (c) 2000-2009 // Joerg Walter, Mathias Koch, Gunter Winkler // // Distributed under the Boost Software License, Version 1.0. (See // accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) // // The authors gratefully acknowledge the support of // GeNeSys mbH & Co. KG in producing this work. // #ifndef _BOOST_UBLAS_FUNCTIONAL_ #define _BOOST_UBLAS_FUNCTIONAL_ #include #include #include #ifdef BOOST_UBLAS_USE_DUFF_DEVICE #include #endif #ifdef BOOST_UBLAS_USE_SIMD #include #else namespace boost { namespace numeric { namespace ublas { namespace raw { }}}} #endif #ifdef BOOST_UBLAS_HAVE_BINDINGS #include #include #include #include #endif #include namespace boost { namespace numeric { namespace ublas { // Scalar functors // Unary template struct scalar_unary_functor { typedef T value_type; typedef typename type_traits::const_reference argument_type; typedef typename type_traits::value_type result_type; }; template struct scalar_identity: public scalar_unary_functor { typedef typename scalar_unary_functor::argument_type argument_type; typedef typename scalar_unary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument_type t) { return t; } }; template struct scalar_negate: public scalar_unary_functor { typedef typename scalar_unary_functor::argument_type argument_type; typedef typename scalar_unary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument_type t) { return - t; } }; template struct scalar_conj: public scalar_unary_functor { typedef typename scalar_unary_functor::value_type value_type; typedef typename scalar_unary_functor::argument_type argument_type; typedef typename scalar_unary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument_type t) { return type_traits::conj (t); } }; // Unary returning real template struct scalar_real_unary_functor { typedef T value_type; typedef typename type_traits::const_reference argument_type; typedef typename type_traits::real_type result_type; }; template struct scalar_real: public scalar_real_unary_functor { typedef typename scalar_real_unary_functor::value_type value_type; typedef typename scalar_real_unary_functor::argument_type argument_type; typedef typename scalar_real_unary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument_type t) { return type_traits::real (t); } }; template struct scalar_imag: public scalar_real_unary_functor { typedef typename scalar_real_unary_functor::value_type value_type; typedef typename scalar_real_unary_functor::argument_type argument_type; typedef typename scalar_real_unary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument_type t) { return type_traits::imag (t); } }; // Binary template struct scalar_binary_functor { typedef typename type_traits::const_reference argument1_type; typedef typename type_traits::const_reference argument2_type; typedef typename promote_traits::promote_type result_type; }; template struct scalar_plus: public scalar_binary_functor { typedef typename scalar_binary_functor::argument1_type argument1_type; typedef typename scalar_binary_functor::argument2_type argument2_type; typedef typename scalar_binary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument1_type t1, argument2_type t2) { return t1 + t2; } }; template struct scalar_minus: public scalar_binary_functor { typedef typename scalar_binary_functor::argument1_type argument1_type; typedef typename scalar_binary_functor::argument2_type argument2_type; typedef typename scalar_binary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument1_type t1, argument2_type t2) { return t1 - t2; } }; template struct scalar_multiplies: public scalar_binary_functor { typedef typename scalar_binary_functor::argument1_type argument1_type; typedef typename scalar_binary_functor::argument2_type argument2_type; typedef typename scalar_binary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument1_type t1, argument2_type t2) { return t1 * t2; } }; template struct scalar_divides: public scalar_binary_functor { typedef typename scalar_binary_functor::argument1_type argument1_type; typedef typename scalar_binary_functor::argument2_type argument2_type; typedef typename scalar_binary_functor::result_type result_type; static BOOST_UBLAS_INLINE result_type apply (argument1_type t1, argument2_type t2) { return t1 / t2; } }; template struct scalar_binary_assign_functor { // ISSUE Remove reference to avoid reference to reference problems typedef typename type_traits::type>::reference argument1_type; typedef typename type_traits::const_reference argument2_type; }; struct assign_tag {}; struct computed_assign_tag {}; template struct scalar_assign: public scalar_binary_assign_functor { typedef typename scalar_binary_assign_functor::argument1_type argument1_type; typedef typename scalar_binary_assign_functor::argument2_type argument2_type; #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) static const bool computed ; #else static const bool computed = false ; #endif static BOOST_UBLAS_INLINE void apply (argument1_type t1, argument2_type t2) { t1 = t2; } template struct rebind { typedef scalar_assign other; }; }; #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) template const bool scalar_assign::computed = false; #endif template struct scalar_plus_assign: public scalar_binary_assign_functor { typedef typename scalar_binary_assign_functor::argument1_type argument1_type; typedef typename scalar_binary_assign_functor::argument2_type argument2_type; #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) static const bool computed ; #else static const bool computed = true ; #endif static BOOST_UBLAS_INLINE void apply (argument1_type t1, argument2_type t2) { t1 += t2; } template struct rebind { typedef scalar_plus_assign other; }; }; #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) template const bool scalar_plus_assign::computed = true; #endif template struct scalar_minus_assign: public scalar_binary_assign_functor { typedef typename scalar_binary_assign_functor::argument1_type argument1_type; typedef typename scalar_binary_assign_functor::argument2_type argument2_type; #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) static const bool computed ; #else static const bool computed = true ; #endif static BOOST_UBLAS_INLINE void apply (argument1_type t1, argument2_type t2) { t1 -= t2; } template struct rebind { typedef scalar_minus_assign other; }; }; #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) template const bool scalar_minus_assign::computed = true; #endif template struct scalar_multiplies_assign: public scalar_binary_assign_functor { typedef typename scalar_binary_assign_functor::argument1_type argument1_type; typedef typename scalar_binary_assign_functor::argument2_type argument2_type; static const bool computed = true; static BOOST_UBLAS_INLINE void apply (argument1_type t1, argument2_type t2) { t1 *= t2; } template struct rebind { typedef scalar_multiplies_assign other; }; }; template struct scalar_divides_assign: public scalar_binary_assign_functor { typedef typename scalar_binary_assign_functor::argument1_type argument1_type; typedef typename scalar_binary_assign_functor::argument2_type argument2_type; static const bool computed ; static BOOST_UBLAS_INLINE void apply (argument1_type t1, argument2_type t2) { t1 /= t2; } template struct rebind { typedef scalar_divides_assign other; }; }; template const bool scalar_divides_assign::computed = true; template struct scalar_binary_swap_functor { typedef typename type_traits::type>::reference argument1_type; typedef typename type_traits::type>::reference argument2_type; }; template struct scalar_swap: public scalar_binary_swap_functor { typedef typename scalar_binary_swap_functor::argument1_type argument1_type; typedef typename scalar_binary_swap_functor::argument2_type argument2_type; static BOOST_UBLAS_INLINE void apply (argument1_type t1, argument2_type t2) { std::swap (t1, t2); } template struct rebind { typedef scalar_swap other; }; }; // Vector functors // Unary returning scalar template struct vector_scalar_unary_functor { typedef typename V::value_type value_type; typedef typename V::value_type result_type; }; template struct vector_sum: public vector_scalar_unary_functor { typedef typename vector_scalar_unary_functor::value_type value_type; typedef typename vector_scalar_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const vector_expression &e) { result_type t = result_type (0); typedef typename E::size_type vector_size_type; vector_size_type size (e ().size ()); for (vector_size_type i = 0; i < size; ++ i) t += e () (i); return t; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (D size, I it) { result_type t = result_type (0); while (-- size >= 0) t += *it, ++ it; return t; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I it, const I &it_end) { result_type t = result_type (0); while (it != it_end) t += *it, ++ it; return t; } }; // Unary returning real scalar template struct vector_scalar_real_unary_functor { typedef typename V::value_type value_type; typedef typename type_traits::real_type real_type; typedef real_type result_type; }; template struct vector_norm_1: public vector_scalar_real_unary_functor { typedef typename vector_scalar_real_unary_functor::value_type value_type; typedef typename vector_scalar_real_unary_functor::real_type real_type; typedef typename vector_scalar_real_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const vector_expression &e) { real_type t = real_type (); typedef typename E::size_type vector_size_type; vector_size_type size (e ().size ()); for (vector_size_type i = 0; i < size; ++ i) { real_type u (type_traits::type_abs (e () (i))); t += u; } return t; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (D size, I it) { real_type t = real_type (); while (-- size >= 0) { real_type u (type_traits::norm_1 (*it)); t += u; ++ it; } return t; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I it, const I &it_end) { real_type t = real_type (); while (it != it_end) { real_type u (type_traits::norm_1 (*it)); t += u; ++ it; } return t; } }; template struct vector_norm_2: public vector_scalar_real_unary_functor { typedef typename vector_scalar_real_unary_functor::value_type value_type; typedef typename vector_scalar_real_unary_functor::real_type real_type; typedef typename vector_scalar_real_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const vector_expression &e) { typedef typename E::size_type vector_size_type; vector_size_type size (e ().size ()); #ifndef BOOST_UBLAS_SCALED_NORM real_type t = real_type (); for (vector_size_type i = 0; i < size; ++ i) { real_type u (type_traits::norm_2 (e () (i))); t += u * u; } return static_cast(type_traits::type_sqrt (t)); #else real_type scale = real_type (); real_type sum_squares (1); for (vector_size_type i = 0; i < size; ++ i) { real_type u (type_traits::norm_2 (e () (i))); if ( real_type () /* zero */ == u ) continue; if (scale < u) { real_type v (scale / u); sum_squares = sum_squares * v * v + real_type (1); scale = u; } else { real_type v (u / scale); sum_squares += v * v; } } return static_cast(scale * type_traits::type_sqrt (sum_squares)); #endif } // Dense case template static BOOST_UBLAS_INLINE result_type apply (D size, I it) { #ifndef BOOST_UBLAS_SCALED_NORM real_type t = real_type (); while (-- size >= 0) { real_type u (type_traits::norm_2 (*it)); t += u * u; ++ it; } return static_cast(type_traits::type_sqrt (t)); #else real_type scale = real_type (); real_type sum_squares (1); while (-- size >= 0) { real_type u (type_traits::norm_2 (*it)); if (scale < u) { real_type v (scale / u); sum_squares = sum_squares * v * v + real_type (1); scale = u; } else { real_type v (u / scale); sum_squares += v * v; } ++ it; } return static_cast(scale * type_traits::type_sqrt (sum_squares)); #endif } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I it, const I &it_end) { #ifndef BOOST_UBLAS_SCALED_NORM real_type t = real_type (); while (it != it_end) { real_type u (type_traits::norm_2 (*it)); t += u * u; ++ it; } return static_cast(type_traits::type_sqrt (t)); #else real_type scale = real_type (); real_type sum_squares (1); while (it != it_end) { real_type u (type_traits::norm_2 (*it)); if (scale < u) { real_type v (scale / u); sum_squares = sum_squares * v * v + real_type (1); scale = u; } else { real_type v (u / scale); sum_squares += v * v; } ++ it; } return static_cast(scale * type_traits::type_sqrt (sum_squares)); #endif } }; template struct vector_norm_2_square : public vector_scalar_real_unary_functor { typedef typename vector_scalar_real_unary_functor::value_type value_type; typedef typename vector_scalar_real_unary_functor::real_type real_type; typedef typename vector_scalar_real_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const vector_expression &e) { real_type t = real_type (); typedef typename E::size_type vector_size_type; vector_size_type size (e ().size ()); for (vector_size_type i = 0; i < size; ++ i) { real_type u (type_traits::norm_2 (e () (i))); t += u * u; } return t; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (D size, I it) { real_type t = real_type (); while (-- size >= 0) { real_type u (type_traits::norm_2 (*it)); t += u * u; ++ it; } return t; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I it, const I &it_end) { real_type t = real_type (); while (it != it_end) { real_type u (type_traits::norm_2 (*it)); t += u * u; ++ it; } return t; } }; template struct vector_norm_inf: public vector_scalar_real_unary_functor { typedef typename vector_scalar_real_unary_functor::value_type value_type; typedef typename vector_scalar_real_unary_functor::real_type real_type; typedef typename vector_scalar_real_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const vector_expression &e) { real_type t = real_type (); typedef typename E::size_type vector_size_type; vector_size_type size (e ().size ()); for (vector_size_type i = 0; i < size; ++ i) { real_type u (type_traits::norm_inf (e () (i))); if (u > t) t = u; } return t; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (D size, I it) { real_type t = real_type (); while (-- size >= 0) { real_type u (type_traits::norm_inf (*it)); if (u > t) t = u; ++ it; } return t; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I it, const I &it_end) { real_type t = real_type (); while (it != it_end) { real_type u (type_traits::norm_inf (*it)); if (u > t) t = u; ++ it; } return t; } }; // Unary returning index template struct vector_scalar_index_unary_functor { typedef typename V::value_type value_type; typedef typename type_traits::real_type real_type; typedef typename V::size_type result_type; }; template struct vector_index_norm_inf: public vector_scalar_index_unary_functor { typedef typename vector_scalar_index_unary_functor::value_type value_type; typedef typename vector_scalar_index_unary_functor::real_type real_type; typedef typename vector_scalar_index_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const vector_expression &e) { // ISSUE For CBLAS compatibility return 0 index in empty case result_type i_norm_inf (0); real_type t = real_type (); typedef typename E::size_type vector_size_type; vector_size_type size (e ().size ()); for (vector_size_type i = 0; i < size; ++ i) { real_type u (type_traits::norm_inf (e () (i))); if (u > t) { i_norm_inf = i; t = u; } } return i_norm_inf; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (D size, I it) { // ISSUE For CBLAS compatibility return 0 index in empty case result_type i_norm_inf (0); real_type t = real_type (); while (-- size >= 0) { real_type u (type_traits::norm_inf (*it)); if (u > t) { i_norm_inf = it.index (); t = u; } ++ it; } return i_norm_inf; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I it, const I &it_end) { // ISSUE For CBLAS compatibility return 0 index in empty case result_type i_norm_inf (0); real_type t = real_type (); while (it != it_end) { real_type u (type_traits::norm_inf (*it)); if (u > t) { i_norm_inf = it.index (); t = u; } ++ it; } return i_norm_inf; } }; // Binary returning scalar template struct vector_scalar_binary_functor { typedef TV value_type; typedef TV result_type; }; template struct vector_inner_prod: public vector_scalar_binary_functor { typedef typename vector_scalar_binary_functor::value_type value_type; typedef typename vector_scalar_binary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const vector_container &c1, const vector_container &c2) { #ifdef BOOST_UBLAS_USE_SIMD using namespace raw; typedef typename C1::size_type vector_size_type; vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ())); const typename V1::value_type *data1 = data_const (c1 ()); const typename V1::value_type *data2 = data_const (c2 ()); vector_size_type s1 = stride (c1 ()); vector_size_type s2 = stride (c2 ()); result_type t = result_type (0); if (s1 == 1 && s2 == 1) { for (vector_size_type i = 0; i < size; ++ i) t += data1 [i] * data2 [i]; } else if (s2 == 1) { for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1) t += data1 [i1] * data2 [i]; } else if (s1 == 1) { for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2) t += data1 [i] * data2 [i2]; } else { for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2) t += data1 [i1] * data2 [i2]; } return t; #elif defined(BOOST_UBLAS_HAVE_BINDINGS) return boost::numeric::bindings::atlas::dot (c1 (), c2 ()); #else return apply (static_cast > (c1), static_cast > (c2)); #endif } template static BOOST_UBLAS_INLINE result_type apply (const vector_expression &e1, const vector_expression &e2) { typedef typename E1::size_type vector_size_type; vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ())); result_type t = result_type (0); #ifndef BOOST_UBLAS_USE_DUFF_DEVICE for (vector_size_type i = 0; i < size; ++ i) t += e1 () (i) * e2 () (i); #else vector_size_type i (0); DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i)); #endif return t; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (D size, I1 it1, I2 it2) { result_type t = result_type (0); #ifndef BOOST_UBLAS_USE_DUFF_DEVICE while (-- size >= 0) t += *it1 * *it2, ++ it1, ++ it2; #else DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); #endif return t; } // Packed case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { result_type t = result_type (0); typedef typename I1::difference_type vector_difference_type; vector_difference_type it1_size (it1_end - it1); vector_difference_type it2_size (it2_end - it2); vector_difference_type diff (0); if (it1_size > 0 && it2_size > 0) diff = it2.index () - it1.index (); if (diff != 0) { vector_difference_type size = (std::min) (diff, it1_size); if (size > 0) { it1 += size; it1_size -= size; diff -= size; } size = (std::min) (- diff, it2_size); if (size > 0) { it2 += size; it2_size -= size; diff += size; } } vector_difference_type size ((std::min) (it1_size, it2_size)); while (-- size >= 0) t += *it1 * *it2, ++ it1, ++ it2; return t; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { result_type t = result_type (0); if (it1 != it1_end && it2 != it2_end) { for (;;) { if (it1.index () == it2.index ()) { t += *it1 * *it2, ++ it1, ++ it2; if (it1 == it1_end || it2 == it2_end) break; } else if (it1.index () < it2.index ()) { increment (it1, it1_end, it2.index () - it1.index ()); if (it1 == it1_end) break; } else if (it1.index () > it2.index ()) { increment (it2, it2_end, it1.index () - it2.index ()); if (it2 == it2_end) break; } } } return t; } }; // Matrix functors // Binary returning vector template struct matrix_vector_binary_functor { typedef typename M1::size_type size_type; typedef typename M1::difference_type difference_type; typedef TV value_type; typedef TV result_type; }; template struct matrix_vector_prod1: public matrix_vector_binary_functor { typedef typename matrix_vector_binary_functor::size_type size_type; typedef typename matrix_vector_binary_functor::difference_type difference_type; typedef typename matrix_vector_binary_functor::value_type value_type; typedef typename matrix_vector_binary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const matrix_container &c1, const vector_container &c2, size_type i) { #ifdef BOOST_UBLAS_USE_SIMD using namespace raw; size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ()); const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ()); const typename M2::value_type *data2 = data_const (c2 ()); size_type s1 = stride2 (c1 ()); size_type s2 = stride (c2 ()); result_type t = result_type (0); if (s1 == 1 && s2 == 1) { for (size_type j = 0; j < size; ++ j) t += data1 [j] * data2 [j]; } else if (s2 == 1) { for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1) t += data1 [j1] * data2 [j]; } else if (s1 == 1) { for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2) t += data1 [j] * data2 [j2]; } else { for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2) t += data1 [j1] * data2 [j2]; } return t; #elif defined(BOOST_UBLAS_HAVE_BINDINGS) return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ()); #else return apply (static_cast > (c1), static_cast > (c2, i)); #endif } template static BOOST_UBLAS_INLINE result_type apply (const matrix_expression &e1, const vector_expression &e2, size_type i) { size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ()); result_type t = result_type (0); #ifndef BOOST_UBLAS_USE_DUFF_DEVICE for (size_type j = 0; j < size; ++ j) t += e1 () (i, j) * e2 () (j); #else size_type j (0); DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j)); #endif return t; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (difference_type size, I1 it1, I2 it2) { result_type t = result_type (0); #ifndef BOOST_UBLAS_USE_DUFF_DEVICE while (-- size >= 0) t += *it1 * *it2, ++ it1, ++ it2; #else DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); #endif return t; } // Packed case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { result_type t = result_type (0); difference_type it1_size (it1_end - it1); difference_type it2_size (it2_end - it2); difference_type diff (0); if (it1_size > 0 && it2_size > 0) diff = it2.index () - it1.index2 (); if (diff != 0) { difference_type size = (std::min) (diff, it1_size); if (size > 0) { it1 += size; it1_size -= size; diff -= size; } size = (std::min) (- diff, it2_size); if (size > 0) { it2 += size; it2_size -= size; diff += size; } } difference_type size ((std::min) (it1_size, it2_size)); while (-- size >= 0) t += *it1 * *it2, ++ it1, ++ it2; return t; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) { result_type t = result_type (0); if (it1 != it1_end && it2 != it2_end) { size_type it1_index = it1.index2 (), it2_index = it2.index (); for (;;) { difference_type compare = it1_index - it2_index; if (compare == 0) { t += *it1 * *it2, ++ it1, ++ it2; if (it1 != it1_end && it2 != it2_end) { it1_index = it1.index2 (); it2_index = it2.index (); } else break; } else if (compare < 0) { increment (it1, it1_end, - compare); if (it1 != it1_end) it1_index = it1.index2 (); else break; } else if (compare > 0) { increment (it2, it2_end, compare); if (it2 != it2_end) it2_index = it2.index (); else break; } } } return t; } // Sparse packed case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */, sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) { result_type t = result_type (0); while (it1 != it1_end) { t += *it1 * it2 () (it1.index2 ()); ++ it1; } return t; } // Packed sparse case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) { result_type t = result_type (0); while (it2 != it2_end) { t += it1 () (it1.index1 (), it2.index ()) * *it2; ++ it2; } return t; } // Another dispatcher template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { typedef typename I1::iterator_category iterator1_category; typedef typename I2::iterator_category iterator2_category; return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ()); } }; template struct matrix_vector_prod2: public matrix_vector_binary_functor { typedef typename matrix_vector_binary_functor::size_type size_type; typedef typename matrix_vector_binary_functor::difference_type difference_type; typedef typename matrix_vector_binary_functor::value_type value_type; typedef typename matrix_vector_binary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const vector_container &c1, const matrix_container &c2, size_type i) { #ifdef BOOST_UBLAS_USE_SIMD using namespace raw; size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ()); const typename M1::value_type *data1 = data_const (c1 ()); const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ()); size_type s1 = stride (c1 ()); size_type s2 = stride1 (c2 ()); result_type t = result_type (0); if (s1 == 1 && s2 == 1) { for (size_type j = 0; j < size; ++ j) t += data1 [j] * data2 [j]; } else if (s2 == 1) { for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1) t += data1 [j1] * data2 [j]; } else if (s1 == 1) { for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2) t += data1 [j] * data2 [j2]; } else { for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2) t += data1 [j1] * data2 [j2]; } return t; #elif defined(BOOST_UBLAS_HAVE_BINDINGS) return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i)); #else return apply (static_cast > (c1), static_cast > (c2, i)); #endif } template static BOOST_UBLAS_INLINE result_type apply (const vector_expression &e1, const matrix_expression &e2, size_type i) { size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ()); result_type t = result_type (0); #ifndef BOOST_UBLAS_USE_DUFF_DEVICE for (size_type j = 0; j < size; ++ j) t += e1 () (j) * e2 () (j, i); #else size_type j (0); DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j)); #endif return t; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (difference_type size, I1 it1, I2 it2) { result_type t = result_type (0); #ifndef BOOST_UBLAS_USE_DUFF_DEVICE while (-- size >= 0) t += *it1 * *it2, ++ it1, ++ it2; #else DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); #endif return t; } // Packed case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { result_type t = result_type (0); difference_type it1_size (it1_end - it1); difference_type it2_size (it2_end - it2); difference_type diff (0); if (it1_size > 0 && it2_size > 0) diff = it2.index1 () - it1.index (); if (diff != 0) { difference_type size = (std::min) (diff, it1_size); if (size > 0) { it1 += size; it1_size -= size; diff -= size; } size = (std::min) (- diff, it2_size); if (size > 0) { it2 += size; it2_size -= size; diff += size; } } difference_type size ((std::min) (it1_size, it2_size)); while (-- size >= 0) t += *it1 * *it2, ++ it1, ++ it2; return t; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) { result_type t = result_type (0); if (it1 != it1_end && it2 != it2_end) { size_type it1_index = it1.index (), it2_index = it2.index1 (); for (;;) { difference_type compare = it1_index - it2_index; if (compare == 0) { t += *it1 * *it2, ++ it1, ++ it2; if (it1 != it1_end && it2 != it2_end) { it1_index = it1.index (); it2_index = it2.index1 (); } else break; } else if (compare < 0) { increment (it1, it1_end, - compare); if (it1 != it1_end) it1_index = it1.index (); else break; } else if (compare > 0) { increment (it2, it2_end, compare); if (it2 != it2_end) it2_index = it2.index1 (); else break; } } } return t; } // Packed sparse case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) { result_type t = result_type (0); while (it2 != it2_end) { t += it1 () (it2.index1 ()) * *it2; ++ it2; } return t; } // Sparse packed case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */, sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) { result_type t = result_type (0); while (it1 != it1_end) { t += *it1 * it2 () (it1.index (), it2.index2 ()); ++ it1; } return t; } // Another dispatcher template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { typedef typename I1::iterator_category iterator1_category; typedef typename I2::iterator_category iterator2_category; return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ()); } }; // Binary returning matrix template struct matrix_matrix_binary_functor { typedef typename M1::size_type size_type; typedef typename M1::difference_type difference_type; typedef TV value_type; typedef TV result_type; }; template struct matrix_matrix_prod: public matrix_matrix_binary_functor { typedef typename matrix_matrix_binary_functor::size_type size_type; typedef typename matrix_matrix_binary_functor::difference_type difference_type; typedef typename matrix_matrix_binary_functor::value_type value_type; typedef typename matrix_matrix_binary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const matrix_container &c1, const matrix_container &c2, size_type i, size_type j) { #ifdef BOOST_UBLAS_USE_SIMD using namespace raw; size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ()); const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ()); const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ()); size_type s1 = stride2 (c1 ()); size_type s2 = stride1 (c2 ()); result_type t = result_type (0); if (s1 == 1 && s2 == 1) { for (size_type k = 0; k < size; ++ k) t += data1 [k] * data2 [k]; } else if (s2 == 1) { for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1) t += data1 [k1] * data2 [k]; } else if (s1 == 1) { for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2) t += data1 [k] * data2 [k2]; } else { for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2) t += data1 [k1] * data2 [k2]; } return t; #elif defined(BOOST_UBLAS_HAVE_BINDINGS) return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j)); #else boost::ignore_unused(j); return apply (static_cast > (c1), static_cast > (c2, i)); #endif } template static BOOST_UBLAS_INLINE result_type apply (const matrix_expression &e1, const matrix_expression &e2, size_type i, size_type j) { size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()); result_type t = result_type (0); #ifndef BOOST_UBLAS_USE_DUFF_DEVICE for (size_type k = 0; k < size; ++ k) t += e1 () (i, k) * e2 () (k, j); #else size_type k (0); DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k)); #endif return t; } // Dense case template static BOOST_UBLAS_INLINE result_type apply (difference_type size, I1 it1, I2 it2) { result_type t = result_type (0); #ifndef BOOST_UBLAS_USE_DUFF_DEVICE while (-- size >= 0) t += *it1 * *it2, ++ it1, ++ it2; #else DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); #endif return t; } // Packed case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) { result_type t = result_type (0); difference_type it1_size (it1_end - it1); difference_type it2_size (it2_end - it2); difference_type diff (0); if (it1_size > 0 && it2_size > 0) diff = it2.index1 () - it1.index2 (); if (diff != 0) { difference_type size = (std::min) (diff, it1_size); if (size > 0) { it1 += size; it1_size -= size; diff -= size; } size = (std::min) (- diff, it2_size); if (size > 0) { it2 += size; it2_size -= size; diff += size; } } difference_type size ((std::min) (it1_size, it2_size)); while (-- size >= 0) t += *it1 * *it2, ++ it1, ++ it2; return t; } // Sparse case template static BOOST_UBLAS_INLINE result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { result_type t = result_type (0); if (it1 != it1_end && it2 != it2_end) { size_type it1_index = it1.index2 (), it2_index = it2.index1 (); for (;;) { difference_type compare = difference_type (it1_index - it2_index); if (compare == 0) { t += *it1 * *it2, ++ it1, ++ it2; if (it1 != it1_end && it2 != it2_end) { it1_index = it1.index2 (); it2_index = it2.index1 (); } else break; } else if (compare < 0) { increment (it1, it1_end, - compare); if (it1 != it1_end) it1_index = it1.index2 (); else break; } else if (compare > 0) { increment (it2, it2_end, compare); if (it2 != it2_end) it2_index = it2.index1 (); else break; } } } return t; } }; // Unary returning scalar norm template struct matrix_scalar_real_unary_functor { typedef typename M::value_type value_type; typedef typename type_traits::real_type real_type; typedef real_type result_type; }; template struct matrix_norm_1: public matrix_scalar_real_unary_functor { typedef typename matrix_scalar_real_unary_functor::value_type value_type; typedef typename matrix_scalar_real_unary_functor::real_type real_type; typedef typename matrix_scalar_real_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const matrix_expression &e) { real_type t = real_type (); typedef typename E::size_type matrix_size_type; matrix_size_type size2 (e ().size2 ()); for (matrix_size_type j = 0; j < size2; ++ j) { real_type u = real_type (); matrix_size_type size1 (e ().size1 ()); for (matrix_size_type i = 0; i < size1; ++ i) { real_type v (type_traits::norm_1 (e () (i, j))); u += v; } if (u > t) t = u; } return t; } }; template struct matrix_norm_frobenius: public matrix_scalar_real_unary_functor { typedef typename matrix_scalar_real_unary_functor::value_type value_type; typedef typename matrix_scalar_real_unary_functor::real_type real_type; typedef typename matrix_scalar_real_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const matrix_expression &e) { real_type t = real_type (); typedef typename E::size_type matrix_size_type; matrix_size_type size1 (e ().size1 ()); for (matrix_size_type i = 0; i < size1; ++ i) { matrix_size_type size2 (e ().size2 ()); for (matrix_size_type j = 0; j < size2; ++ j) { real_type u (type_traits::norm_2 (e () (i, j))); t += u * u; } } return type_traits::type_sqrt (t); } }; template struct matrix_norm_inf: public matrix_scalar_real_unary_functor { typedef typename matrix_scalar_real_unary_functor::value_type value_type; typedef typename matrix_scalar_real_unary_functor::real_type real_type; typedef typename matrix_scalar_real_unary_functor::result_type result_type; template static BOOST_UBLAS_INLINE result_type apply (const matrix_expression &e) { real_type t = real_type (); typedef typename E::size_type matrix_size_type; matrix_size_type size1 (e ().size1 ()); for (matrix_size_type i = 0; i < size1; ++ i) { real_type u = real_type (); matrix_size_type size2 (e ().size2 ()); for (matrix_size_type j = 0; j < size2; ++ j) { real_type v (type_traits::norm_inf (e () (i, j))); u += v; } if (u > t) t = u; } return t; } }; // forward declaration template struct basic_column_major; // This functor defines storage layout and it's properties // matrix (i,j) -> storage [i * size_i + j] template struct basic_row_major { typedef Z size_type; typedef D difference_type; typedef row_major_tag orientation_category; typedef basic_column_major transposed_layout; static BOOST_UBLAS_INLINE size_type storage_size (size_type size_i, size_type size_j) { // Guard against size_type overflow BOOST_UBLAS_CHECK (size_j == 0 || size_i <= (std::numeric_limits::max) () / size_j, bad_size ()); return size_i * size_j; } // Indexing conversion to storage element static BOOST_UBLAS_INLINE size_type element (size_type i, size_type size_i, size_type j, size_type size_j) { BOOST_UBLAS_CHECK (i < size_i, bad_index ()); BOOST_UBLAS_CHECK (j < size_j, bad_index ()); detail::ignore_unused_variable_warning(size_i); // Guard against size_type overflow BOOST_UBLAS_CHECK (i <= ((std::numeric_limits::max) () - j) / size_j, bad_index ()); return i * size_j + j; } static BOOST_UBLAS_INLINE size_type address (size_type i, size_type size_i, size_type j, size_type size_j) { BOOST_UBLAS_CHECK (i <= size_i, bad_index ()); BOOST_UBLAS_CHECK (j <= size_j, bad_index ()); // Guard against size_type overflow - address may be size_j past end of storage BOOST_UBLAS_CHECK (size_j == 0 || i <= ((std::numeric_limits::max) () - j) / size_j, bad_index ()); detail::ignore_unused_variable_warning(size_i); return i * size_j + j; } // Storage element to index conversion static BOOST_UBLAS_INLINE difference_type distance_i (difference_type k, size_type /* size_i */, size_type size_j) { return size_j != 0 ? k / size_j : 0; } static BOOST_UBLAS_INLINE difference_type distance_j (difference_type k, size_type /* size_i */, size_type /* size_j */) { return k; } static BOOST_UBLAS_INLINE size_type index_i (difference_type k, size_type /* size_i */, size_type size_j) { return size_j != 0 ? k / size_j : 0; } static BOOST_UBLAS_INLINE size_type index_j (difference_type k, size_type /* size_i */, size_type size_j) { return size_j != 0 ? k % size_j : 0; } static BOOST_UBLAS_INLINE bool fast_i () { return false; } static BOOST_UBLAS_INLINE bool fast_j () { return true; } // Iterating storage elements template static BOOST_UBLAS_INLINE void increment_i (I &it, size_type /* size_i */, size_type size_j) { it += size_j; } template static BOOST_UBLAS_INLINE void increment_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) { it += n * size_j; } template static BOOST_UBLAS_INLINE void decrement_i (I &it, size_type /* size_i */, size_type size_j) { it -= size_j; } template static BOOST_UBLAS_INLINE void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) { it -= n * size_j; } template static BOOST_UBLAS_INLINE void increment_j (I &it, size_type /* size_i */, size_type /* size_j */) { ++ it; } template static BOOST_UBLAS_INLINE void increment_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { it += n; } template static BOOST_UBLAS_INLINE void decrement_j (I &it, size_type /* size_i */, size_type /* size_j */) { -- it; } template static BOOST_UBLAS_INLINE void decrement_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { it -= n; } // Triangular access static BOOST_UBLAS_INLINE size_type triangular_size (size_type size_i, size_type size_j) { size_type size = (std::max) (size_i, size_j); // Guard against size_type overflow - simplified BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits::max) () / size /* +1/2 */, bad_size ()); return ((size + 1) * size) / 2; } static BOOST_UBLAS_INLINE size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) { BOOST_UBLAS_CHECK (i < size_i, bad_index ()); BOOST_UBLAS_CHECK (j < size_j, bad_index ()); BOOST_UBLAS_CHECK (i >= j, bad_index ()); detail::ignore_unused_variable_warning(size_i); detail::ignore_unused_variable_warning(size_j); // FIXME size_type overflow // sigma_i (i + 1) = (i + 1) * i / 2 // i = 0 1 2 3, sigma = 0 1 3 6 return ((i + 1) * i) / 2 + j; } static BOOST_UBLAS_INLINE size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) { BOOST_UBLAS_CHECK (i < size_i, bad_index ()); BOOST_UBLAS_CHECK (j < size_j, bad_index ()); BOOST_UBLAS_CHECK (i <= j, bad_index ()); // FIXME size_type overflow // sigma_i (size - i) = size * i - i * (i - 1) / 2 // i = 0 1 2 3, sigma = 0 4 7 9 return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i; } // Major and minor indices static BOOST_UBLAS_INLINE size_type index_M (size_type index1, size_type /* index2 */) { return index1; } static BOOST_UBLAS_INLINE size_type index_m (size_type /* index1 */, size_type index2) { return index2; } static BOOST_UBLAS_INLINE size_type size_M (size_type size_i, size_type /* size_j */) { return size_i; } static BOOST_UBLAS_INLINE size_type size_m (size_type /* size_i */, size_type size_j) { return size_j; } }; // This functor defines storage layout and it's properties // matrix (i,j) -> storage [i + j * size_i] template struct basic_column_major { typedef Z size_type; typedef D difference_type; typedef column_major_tag orientation_category; typedef basic_row_major transposed_layout; static BOOST_UBLAS_INLINE size_type storage_size (size_type size_i, size_type size_j) { // Guard against size_type overflow BOOST_UBLAS_CHECK (size_i == 0 || size_j <= (std::numeric_limits::max) () / size_i, bad_size ()); return size_i * size_j; } // Indexing conversion to storage element static BOOST_UBLAS_INLINE size_type element (size_type i, size_type size_i, size_type j, size_type size_j) { BOOST_UBLAS_CHECK (i < size_i, bad_index ()); BOOST_UBLAS_CHECK (j < size_j, bad_index ()); detail::ignore_unused_variable_warning(size_j); // Guard against size_type overflow BOOST_UBLAS_CHECK (j <= ((std::numeric_limits::max) () - i) / size_i, bad_index ()); return i + j * size_i; } static BOOST_UBLAS_INLINE size_type address (size_type i, size_type size_i, size_type j, size_type size_j) { BOOST_UBLAS_CHECK (i <= size_i, bad_index ()); BOOST_UBLAS_CHECK (j <= size_j, bad_index ()); detail::ignore_unused_variable_warning(size_j); // Guard against size_type overflow - address may be size_i past end of storage BOOST_UBLAS_CHECK (size_i == 0 || j <= ((std::numeric_limits::max) () - i) / size_i, bad_index ()); return i + j * size_i; } // Storage element to index conversion static BOOST_UBLAS_INLINE difference_type distance_i (difference_type k, size_type /* size_i */, size_type /* size_j */) { return k; } static BOOST_UBLAS_INLINE difference_type distance_j (difference_type k, size_type size_i, size_type /* size_j */) { return size_i != 0 ? k / size_i : 0; } static BOOST_UBLAS_INLINE size_type index_i (difference_type k, size_type size_i, size_type /* size_j */) { return size_i != 0 ? k % size_i : 0; } static BOOST_UBLAS_INLINE size_type index_j (difference_type k, size_type size_i, size_type /* size_j */) { return size_i != 0 ? k / size_i : 0; } static BOOST_UBLAS_INLINE bool fast_i () { return true; } static BOOST_UBLAS_INLINE bool fast_j () { return false; } // Iterating template static BOOST_UBLAS_INLINE void increment_i (I &it, size_type /* size_i */, size_type /* size_j */) { ++ it; } template static BOOST_UBLAS_INLINE void increment_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { it += n; } template static BOOST_UBLAS_INLINE void decrement_i (I &it, size_type /* size_i */, size_type /* size_j */) { -- it; } template static BOOST_UBLAS_INLINE void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { it -= n; } template static BOOST_UBLAS_INLINE void increment_j (I &it, size_type size_i, size_type /* size_j */) { it += size_i; } template static BOOST_UBLAS_INLINE void increment_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) { it += n * size_i; } template static BOOST_UBLAS_INLINE void decrement_j (I &it, size_type size_i, size_type /* size_j */) { it -= size_i; } template static BOOST_UBLAS_INLINE void decrement_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) { it -= n* size_i; } // Triangular access static BOOST_UBLAS_INLINE size_type triangular_size (size_type size_i, size_type size_j) { size_type size = (std::max) (size_i, size_j); // Guard against size_type overflow - simplified BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits::max) () / size /* +1/2 */, bad_size ()); return ((size + 1) * size) / 2; } static BOOST_UBLAS_INLINE size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) { BOOST_UBLAS_CHECK (i < size_i, bad_index ()); BOOST_UBLAS_CHECK (j < size_j, bad_index ()); BOOST_UBLAS_CHECK (i >= j, bad_index ()); // FIXME size_type overflow // sigma_j (size - j) = size * j - j * (j - 1) / 2 // j = 0 1 2 3, sigma = 0 4 7 9 return i - j + (j * (2 * (std::max) (size_i, size_j) - j + 1)) / 2; } static BOOST_UBLAS_INLINE size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) { BOOST_UBLAS_CHECK (i < size_i, bad_index ()); BOOST_UBLAS_CHECK (j < size_j, bad_index ()); BOOST_UBLAS_CHECK (i <= j, bad_index ()); boost::ignore_unused(size_i, size_j); // FIXME size_type overflow // sigma_j (j + 1) = (j + 1) * j / 2 // j = 0 1 2 3, sigma = 0 1 3 6 return i + ((j + 1) * j) / 2; } // Major and minor indices static BOOST_UBLAS_INLINE size_type index_M (size_type /* index1 */, size_type index2) { return index2; } static BOOST_UBLAS_INLINE size_type index_m (size_type index1, size_type /* index2 */) { return index1; } static BOOST_UBLAS_INLINE size_type size_M (size_type /* size_i */, size_type size_j) { return size_j; } static BOOST_UBLAS_INLINE size_type size_m (size_type size_i, size_type /* size_j */) { return size_i; } }; template struct basic_full { typedef Z size_type; template static BOOST_UBLAS_INLINE size_type packed_size (L, size_type size_i, size_type size_j) { return L::storage_size (size_i, size_j); } static BOOST_UBLAS_INLINE bool zero (size_type /* i */, size_type /* j */) { return false; } static BOOST_UBLAS_INLINE bool one (size_type /* i */, size_type /* j */) { return false; } static BOOST_UBLAS_INLINE bool other (size_type /* i */, size_type /* j */) { return true; } // FIXME: this should not be used at all static BOOST_UBLAS_INLINE size_type restrict1 (size_type i, size_type /* j */) { return i; } static BOOST_UBLAS_INLINE size_type restrict2 (size_type /* i */, size_type j) { return j; } static BOOST_UBLAS_INLINE size_type mutable_restrict1 (size_type i, size_type /* j */) { return i; } static BOOST_UBLAS_INLINE size_type mutable_restrict2 (size_type /* i */, size_type j) { return j; } }; namespace detail { template < class L > struct transposed_structure { typedef typename L::size_type size_type; template static BOOST_UBLAS_INLINE size_type packed_size (LAYOUT l, size_type size_i, size_type size_j) { return L::packed_size(l, size_j, size_i); } static BOOST_UBLAS_INLINE bool zero (size_type i, size_type j) { return L::zero(j, i); } static BOOST_UBLAS_INLINE bool one (size_type i, size_type j) { return L::one(j, i); } static BOOST_UBLAS_INLINE bool other (size_type i, size_type j) { return L::other(j, i); } template static BOOST_UBLAS_INLINE size_type element (LAYOUT /* l */, size_type i, size_type size_i, size_type j, size_type size_j) { return L::element(typename LAYOUT::transposed_layout(), j, size_j, i, size_i); } static BOOST_UBLAS_INLINE size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { return L::restrict2(j, i, size2, size1); } static BOOST_UBLAS_INLINE size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) { return L::restrict1(j, i, size2, size1); } static BOOST_UBLAS_INLINE size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type size2) { return L::mutable_restrict2(j, i, size2, size1); } static BOOST_UBLAS_INLINE size_type mutable_restrict2 (size_type i, size_type j, size_type size1, size_type size2) { return L::mutable_restrict1(j, i, size2, size1); } static BOOST_UBLAS_INLINE size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { return L::global_restrict2(index2, size2, index1, size1); } static BOOST_UBLAS_INLINE size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { return L::global_restrict1(index2, size2, index1, size1); } static BOOST_UBLAS_INLINE size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { return L::global_mutable_restrict2(index2, size2, index1, size1); } static BOOST_UBLAS_INLINE size_type global_mutable_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { return L::global_mutable_restrict1(index2, size2, index1, size1); } }; } template struct basic_lower { typedef Z size_type; typedef lower_tag triangular_type; template static BOOST_UBLAS_INLINE size_type packed_size (L, size_type size_i, size_type size_j) { return L::triangular_size (size_i, size_j); } static BOOST_UBLAS_INLINE bool zero (size_type i, size_type j) { return j > i; } static BOOST_UBLAS_INLINE bool one (size_type /* i */, size_type /* j */) { return false; } static BOOST_UBLAS_INLINE bool other (size_type i, size_type j) { return j <= i; } template static BOOST_UBLAS_INLINE size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { return L::lower_element (i, size_i, j, size_j); } // return nearest valid index in column j static BOOST_UBLAS_INLINE size_type restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { return (std::max)(j, (std::min) (size1, i)); } // return nearest valid index in row i static BOOST_UBLAS_INLINE size_type restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { return (std::max)(size_type(0), (std::min) (i+1, j)); } // return nearest valid mutable index in column j static BOOST_UBLAS_INLINE size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { return (std::max)(j, (std::min) (size1, i)); } // return nearest valid mutable index in row i static BOOST_UBLAS_INLINE size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { return (std::max)(size_type(0), (std::min) (i+1, j)); } // return an index between the first and (1+last) filled row static BOOST_UBLAS_INLINE size_type global_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { return (std::max)(size_type(0), (std::min)(size1, index1) ); } // return an index between the first and (1+last) filled column static BOOST_UBLAS_INLINE size_type global_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { return (std::max)(size_type(0), (std::min)(size2, index2) ); } // return an index between the first and (1+last) filled mutable row static BOOST_UBLAS_INLINE size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { return (std::max)(size_type(0), (std::min)(size1, index1) ); } // return an index between the first and (1+last) filled mutable column static BOOST_UBLAS_INLINE size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { return (std::max)(size_type(0), (std::min)(size2, index2) ); } }; // the first row only contains a single 1. Thus it is not stored. template struct basic_unit_lower : public basic_lower { typedef Z size_type; typedef unit_lower_tag triangular_type; template static BOOST_UBLAS_INLINE size_type packed_size (L, size_type size_i, size_type size_j) { // Zero size strict triangles are bad at this point BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ()); return L::triangular_size (size_i - 1, size_j - 1); } static BOOST_UBLAS_INLINE bool one (size_type i, size_type j) { return j == i; } static BOOST_UBLAS_INLINE bool other (size_type i, size_type j) { return j < i; } template static BOOST_UBLAS_INLINE size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { // Zero size strict triangles are bad at this point BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ()); return L::lower_element (i-1, size_i - 1, j, size_j - 1); } static BOOST_UBLAS_INLINE size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { return (std::max)(j+1, (std::min) (size1, i)); } static BOOST_UBLAS_INLINE size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { return (std::max)(size_type(0), (std::min) (i, j)); } // return an index between the first and (1+last) filled mutable row static BOOST_UBLAS_INLINE size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { return (std::max)(size_type(1), (std::min)(size1, index1) ); } // return an index between the first and (1+last) filled mutable column static BOOST_UBLAS_INLINE size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { BOOST_UBLAS_CHECK( size2 >= 1 , external_logic() ); return (std::max)(size_type(0), (std::min)(size2-1, index2) ); } }; // the first row only contains no element. Thus it is not stored. template struct basic_strict_lower : public basic_unit_lower { typedef Z size_type; typedef strict_lower_tag triangular_type; template static BOOST_UBLAS_INLINE size_type packed_size (L, size_type size_i, size_type size_j) { // Zero size strict triangles are bad at this point BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ()); return L::triangular_size (size_i - 1, size_j - 1); } static BOOST_UBLAS_INLINE bool zero (size_type i, size_type j) { return j >= i; } static BOOST_UBLAS_INLINE bool one (size_type /*i*/, size_type /*j*/) { return false; } static BOOST_UBLAS_INLINE bool other (size_type i, size_type j) { return j < i; } template static BOOST_UBLAS_INLINE size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { // Zero size strict triangles are bad at this point BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ()); return L::lower_element (i-1, size_i - 1, j, size_j - 1); } static BOOST_UBLAS_INLINE size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { return basic_unit_lower::mutable_restrict1(i, j, size1, size2); } static BOOST_UBLAS_INLINE size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) { return basic_unit_lower::mutable_restrict2(i, j, size1, size2); } // return an index between the first and (1+last) filled row static BOOST_UBLAS_INLINE size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { return basic_unit_lower::global_mutable_restrict1(index1, size1, index2, size2); } // return an index between the first and (1+last) filled column static BOOST_UBLAS_INLINE size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { return basic_unit_lower::global_mutable_restrict2(index1, size1, index2, size2); } }; template struct basic_upper : public detail::transposed_structure > { typedef upper_tag triangular_type; }; template struct basic_unit_upper : public detail::transposed_structure > { typedef unit_upper_tag triangular_type; }; template struct basic_strict_upper : public detail::transposed_structure > { typedef strict_upper_tag triangular_type; }; }}} #endif