functional.hpp

boost库提供标准的C++ API 配合dev c++使用,功能更加强大

            size_type s2 = stride (e2 ());
            result_type t = result_type ();
            if (s1 == 1 && s2 == 1) {
                for (size_type i = 0; i < size; ++ i)
                    t += data1 [i] * data2 [i];
            } else if (s2 == 1) {
                for (size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
                    t += data1 [i1] * data2 [i];
            } else if (s1 == 1) {
                for (size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
                    t += data1 [i] * data2 [i2];
            } else {
                for (size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
                    t += data1 [i1] * data2 [i2];
            }
            return t;
#else
            return boost::numeric::bindings::atlas::dot (e1 (), e2 ());
#endif
        }
        template<class E1, class E2>
        BOOST_UBLAS_INLINE
        result_type operator () (const vector_expression<E1> &e1,
                                 const vector_expression<E2> &e2,
                                 abstract_tag) const {
            size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
            result_type t = result_type ();
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            for (size_type i = 0; i < size; ++ i)
                t += e1 () (i) * e2 () (i);
#else
            size_type i (0);
            DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));
#endif
            return t;
        }
        template<class E1, class E2>
        BOOST_UBLAS_INLINE
        result_type operator () (const vector_expression<E1> &e1,
                                 const vector_expression<E2> &e2) const {
#ifdef BOOST_UBLAS_USE_SIMD
            typedef typename boost::mpl::if_c<
                boost::mpl::and_<boost::is_same<typename E1::simd_category, concrete_tag>,
                                 boost::is_same<typename E2::simd_category, concrete_tag> >::value,
                    concrete_tag,
                    abstract_tag>::type simd_category;
#else
            typedef abstract_tag simd_category;
#endif
            return operator () (e1, e2, simd_category ());
        }
        // Dense case
        template<class I1, class I2>
        BOOST_UBLAS_INLINE
        result_type operator () (difference_type size, I1 it1, I2 it2) const {
            result_type t = result_type ();
#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>
        BOOST_UBLAS_INLINE
        result_type operator () (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) const {
            result_type t = result_type ();
            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.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>
        BOOST_UBLAS_INLINE
        result_type operator () (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) const {
            result_type t = result_type ();
            while (it1 != it1_end && it2 != it2_end) {
                difference_type compare = it1.index () - it2.index ();
                if (compare < 0) 
                    ++ it1;
                else if (compare == 0) 
                    t += *it1 * *it2, ++ it1, ++ it2;
                else if (compare > 0)
                    ++ it2;
            }
            return t; 
        }
    };

    // Matrix functors

    // Binary returning vector
    template<class T1, class T2, class TR>
    struct matrix_vector_binary_functor {
        typedef std::size_t size_type;
        typedef std::ptrdiff_t difference_type;
        typedef TR value_type;
        typedef TR result_type;
    };

    template<class T1, class T2, class TR>
    struct matrix_vector_prod1:
        public matrix_vector_binary_functor<T1, T2, TR> {
        typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;
        typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;
        typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;
        typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;

        template<class E1, class E2>
        BOOST_UBLAS_INLINE
        result_type operator () (const matrix_expression<E1> &e1,
                                 const vector_expression<E2> &e2,
                                 size_type i, concrete_tag) const {
#ifndef BOOST_UBLAS_HAVE_BINDINGS
            size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
            const T1 *data1 = data_const (e1 ()) + i * stride1 (e1 ());
            const T2 *data2 = data_const (e2 ());
            size_type s1 = stride2 (e1 ());
            size_type s2 = stride (e2 ());
            result_type t = result_type ();
            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;
#else
            return boost::numeric::bindings::atlas::dot (e1 ().row (i), e2 ());
#endif
        }
        template<class E1, class E2>
        BOOST_UBLAS_INLINE
        result_type operator () (const matrix_expression<E1> &e1,
                                 const vector_expression<E2> &e2,
                                 size_type i, abstract_tag) const {
            size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
            result_type t = result_type ();
#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;
        }
        template<class E1, class E2>
        BOOST_UBLAS_INLINE
        result_type operator () (const matrix_expression<E1> &e1,
                                 const vector_expression<E2> &e2,
                                 size_type i) const {
#ifdef BOOST_UBLAS_USE_SIMD
            typedef typename boost::mpl::if_c<
                boost::mpl::and_<boost::is_same<typename E1::simd_category, concrete_tag>,
                                 boost::is_same<typename E2::simd_category, concrete_tag> >::value,
                    concrete_tag,
                    abstract_tag>::type simd_category;
#else
            typedef abstract_tag simd_category;
#endif
            return operator () (e1, e2, i, simd_category ());
        }
        // Dense case
        template<class I1, class I2>
        BOOST_UBLAS_INLINE
        result_type operator () (difference_type size, I1 it1, I2 it2) const {
            result_type t = result_type ();
#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>
        BOOST_UBLAS_INLINE
        result_type operator () (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) const {
            result_type t = result_type ();
            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>
        BOOST_UBLAS_INLINE
        result_type operator () (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
                                 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) const {
            result_type t = result_type ();
            while (it1 != it1_end && it2 != it2_end) {
                difference_type compare = it1.index2 () - it2.index ();
                if (compare < 0)
                    ++ it1;
                else if (compare == 0)
                    t += *it1 * *it2, ++ it1, ++ it2;
                else if (compare > 0)
                    ++ it2;
            }
            return t;
        }
        // Sparse packed case
        template<class I1, class I2>
        BOOST_UBLAS_INLINE
        result_type operator () (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
                                 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) const {
            result_type t = result_type ();
            while (it1 != it1_end) {
                t += *it1 * it2 () (it1.index2 ());
                ++ it1;
            }
            return t;
        }
        // Packed sparse case
        template<class I1, class I2>
        BOOST_UBLAS_INLINE
        result_type operator () (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
                                 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) const {
            result_type t = result_type ();
            while (it2 != it2_end) {
                t += it1 () (it1.index1 (), it2.index ()) * *it2;
                ++ it2;
            }
            return t;
        }
        // Another dispatcher
        template<class I1, class I2>
        BOOST_UBLAS_INLINE
        result_type operator () (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
                                 sparse_bidirectional_iterator_tag) const {
            typedef typename I1::iterator_category iterator1_category;
            typedef typename I2::iterator_category iterator2_category;
            return operator () (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
        }
    };

    template<class T1, class T2, class TR>
    struct matrix_vector_prod2:
        public matrix_vector_binary_functor<T1, T2, TR> {
        typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;
        typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;
        typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;
        typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;

        template<class E1, class E2>
        BOOST_UBLAS_INLINE
        result_type operator () (const vector_expression<E1> &e1,
                                 const matrix_expression<E2> &e2,
                                 size_type i, concrete_tag) const {
#ifndef BOOST_UBLAS_HAVE_BINDINGS
            size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
            const T1 *data1 = data_const (e1 ());
            const T2 *data2 = data_const (e2 ()) + i * stride2 (e2 ());
            size_type s1 = stride (e1 ());
            size_type s2 = stride1 (e2 ());
            result_type t = result_type ();
            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;
#else
            return boost::numeric::bindings::atlas::dot (e1 (), e2 ().column (i));
#endif
        }
        template<class E1, class E2>
        BOOST_UBLAS_INLINE
        result_type operator () (const vector_expression<E1> &e1,