source: NonGTP/Boost/boost/numeric/ublas/vector_expression.hpp @ 857

Revision 857, 61.2 KB checked in by igarcia, 19 years ago (diff)
Line 
1//
2//  Copyright (c) 2000-2002
3//  Joerg Walter, Mathias Koch
4//
5//  Permission to use, copy, modify, distribute and sell this software
6//  and its documentation for any purpose is hereby granted without fee,
7//  provided that the above copyright notice appear in all copies and
8//  that both that copyright notice and this permission notice appear
9//  in supporting documentation.  The authors make no representations
10//  about the suitability of this software for any purpose.
11//  It is provided "as is" without express or implied warranty.
12//
13//  The authors gratefully acknowledge the support of
14//  GeNeSys mbH & Co. KG in producing this work.
15//
16
17#ifndef _BOOST_UBLAS_VECTOR_EXPRESSION_
18#define _BOOST_UBLAS_VECTOR_EXPRESSION_
19
20#include <boost/numeric/ublas/expression_types.hpp>
21
22
23// Expression templates based on ideas of Todd Veldhuizen and Geoffrey Furnish
24// Iterators based on ideas of Jeremy Siek
25//
26// Classes that model the Vector Expression concept
27
28namespace boost { namespace numeric { namespace ublas {
29
30    template<class E>
31    class vector_reference:
32        public vector_expression<vector_reference<E> > {
33
34        typedef vector_reference<E> self_type;
35    public:
36#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
37        using vector_expression<vector_reference<E> >::operator ();
38#endif
39        typedef typename E::size_type size_type;
40        typedef typename E::difference_type difference_type;
41        typedef typename E::value_type value_type;
42        typedef typename E::const_reference const_reference;
43        typedef typename boost::mpl::if_<boost::is_const<E>,
44                                          typename E::const_reference,
45                                          typename E::reference>::type reference;
46        typedef E refered_type;
47        typedef const self_type const_closure_type;
48        typedef const_closure_type closure_type;
49        typedef typename E::storage_category storage_category;
50
51        // Construction and destruction
52        BOOST_UBLAS_INLINE
53        explicit vector_reference (refered_type &e):
54            e_ (e) {}
55
56        // Accessors
57        BOOST_UBLAS_INLINE
58        size_type size () const {
59            return expression ().size ();
60        }
61
62    public:
63        // Expression accessors - const correct
64        BOOST_UBLAS_INLINE
65        const refered_type &expression () const {
66            return e_;
67        }
68        BOOST_UBLAS_INLINE
69        refered_type &expression () {
70            return e_;
71        }
72
73    public:
74        // Element access
75#ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
76        BOOST_UBLAS_INLINE
77        const_reference operator () (size_type i) const {
78            return expression () (i);
79        }
80        BOOST_UBLAS_INLINE
81        reference operator () (size_type i) {
82            return expression () (i);
83        }
84
85        BOOST_UBLAS_INLINE
86        const_reference operator [] (size_type i) const {
87            return expression () [i];
88        }
89        BOOST_UBLAS_INLINE
90        reference operator [] (size_type i) {
91            return expression () [i];
92        }
93#else
94        BOOST_UBLAS_INLINE
95        reference operator () (size_type i) const {
96            return expression () (i);
97        }
98
99        BOOST_UBLAS_INLINE
100        reference operator [] (size_type i) const {
101            return expression () [i];
102        }
103#endif
104
105        // Assignment
106        BOOST_UBLAS_INLINE
107        vector_reference &operator = (const vector_reference &v) {
108            expression ().operator = (v);
109            return *this;
110        }
111        template<class AE>
112        BOOST_UBLAS_INLINE
113        vector_reference &operator = (const vector_expression<AE> &ae) {
114            expression ().operator = (ae);
115            return *this;
116        }
117        template<class AE>
118        BOOST_UBLAS_INLINE
119        vector_reference &assign (const vector_expression<AE> &ae) {
120            expression ().assign (ae);
121            return *this;
122        }
123        template<class AE>
124        BOOST_UBLAS_INLINE
125        vector_reference &operator += (const vector_expression<AE> &ae) {
126            expression ().operator += (ae);
127            return *this;
128        }
129        template<class AE>
130        BOOST_UBLAS_INLINE
131        vector_reference &plus_assign (const vector_expression<AE> &ae) {
132            expression ().plus_assign (ae);
133            return *this;
134        }
135        template<class AE>
136        BOOST_UBLAS_INLINE
137        vector_reference &operator -= (const vector_expression<AE> &ae) {
138            expression ().operator -= (ae);
139            return *this;
140        }
141        template<class AE>
142        BOOST_UBLAS_INLINE
143        vector_reference &minus_assign (const vector_expression<AE> &ae) {
144            expression ().minus_assign (ae);
145            return *this;
146        }
147        template<class AT>
148        BOOST_UBLAS_INLINE
149        vector_reference &operator *= (const AT &at) {
150            expression ().operator *= (at);
151            return *this;
152        }
153        template<class AT>
154        BOOST_UBLAS_INLINE
155        vector_reference &operator /= (const AT &at) {
156            expression ().operator /= (at);
157            return *this;
158        }
159
160        // Swapping
161        BOOST_UBLAS_INLINE
162        void swap (vector_reference &v) {
163            expression ().swap (v.expression ());
164        }
165
166        // Closure comparison
167        BOOST_UBLAS_INLINE
168        bool same_closure (const vector_reference &vr) const {
169            return &(*this).e_ == &vr.e_;
170        }
171
172        // Iterator types
173        typedef typename E::const_iterator const_iterator;
174        typedef typename boost::mpl::if_<boost::is_const<E>,
175                                          typename E::const_iterator,
176                                          typename E::iterator>::type iterator;
177
178        // Element lookup
179        BOOST_UBLAS_INLINE
180        const_iterator find (size_type i) const {
181            return expression ().find (i);
182        }
183        BOOST_UBLAS_INLINE
184        iterator find (size_type i) {
185            return expression ().find (i);
186        }
187
188        // Iterator is the iterator of the referenced expression.
189
190        BOOST_UBLAS_INLINE
191        const_iterator begin () const {
192            return expression ().begin ();
193        }
194        BOOST_UBLAS_INLINE
195        const_iterator end () const {
196            return expression ().end ();
197        }
198
199        BOOST_UBLAS_INLINE
200        iterator begin () {
201            return expression ().begin ();
202        }
203        BOOST_UBLAS_INLINE
204        iterator end () {
205            return expression ().end ();
206        }
207
208        // Reverse iterator
209        typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
210        typedef reverse_iterator_base<iterator> reverse_iterator;
211
212        BOOST_UBLAS_INLINE
213        const_reverse_iterator rbegin () const {
214            return const_reverse_iterator (end ());
215        }
216        BOOST_UBLAS_INLINE
217        const_reverse_iterator rend () const {
218            return const_reverse_iterator (begin ());
219        }
220        BOOST_UBLAS_INLINE
221        reverse_iterator rbegin () {
222            return reverse_iterator (end ());
223        }
224        BOOST_UBLAS_INLINE
225        reverse_iterator rend () {
226            return reverse_iterator (begin ());
227        }
228
229    private:
230        refered_type &e_;
231    };
232
233
234    template<class E, class F>
235    class vector_unary:
236        public vector_expression<vector_unary<E, F> > {
237
238        typedef F functor_type;
239        typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<typename E::value_type> >,
240                                          E,
241                                          const E>::type expression_type;
242        typedef typename boost::mpl::if_<boost::is_const<expression_type>,
243                                          typename E::const_closure_type,
244                                          typename E::closure_type>::type expression_closure_type;
245        typedef vector_unary<E, F> self_type;
246    public:
247#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
248        using vector_expression<vector_unary<E, F> >::operator ();
249#endif
250        typedef typename E::size_type size_type;
251        typedef typename E::difference_type difference_type;
252        typedef typename F::result_type value_type;
253        typedef value_type const_reference;
254        typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<value_type> >,
255                                          typename E::reference,
256                                          value_type>::type reference;
257        typedef const self_type const_closure_type;
258        typedef self_type closure_type;
259        typedef unknown_storage_tag storage_category;
260
261        // Construction and destruction
262        BOOST_UBLAS_INLINE
263        // May be used as mutable expression.
264        explicit vector_unary (expression_type &e):
265            e_ (e) {}
266
267        // Accessors
268        BOOST_UBLAS_INLINE
269        size_type size () const {
270            return e_.size ();
271        }
272
273    public:
274        // Expression accessors
275        BOOST_UBLAS_INLINE
276        const expression_closure_type &expression () const {
277            return e_;
278        }
279
280    public:
281        // Element access
282        BOOST_UBLAS_INLINE
283        const_reference operator () (size_type i) const {
284            return functor_type::apply (e_ (i));
285        }
286        BOOST_UBLAS_INLINE
287        reference operator () (size_type i) {
288            BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
289            return e_ (i);
290        }
291
292        BOOST_UBLAS_INLINE
293        const_reference operator [] (size_type i) const {
294            return functor_type::apply (e_ [i]);
295        }
296        BOOST_UBLAS_INLINE
297        reference operator [] (size_type i) {
298            BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
299            return e_ [i];
300        }
301
302        // Closure comparison
303        BOOST_UBLAS_INLINE
304        bool same_closure (const vector_unary &vu) const {
305            return (*this).expression ().same_closure (vu.expression ());
306        }
307
308        // Iterator types
309    private:
310        typedef typename E::const_iterator const_subiterator_type;
311        typedef const value_type *const_pointer;
312
313    public:
314#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
315        typedef indexed_const_iterator<const_closure_type, typename const_subiterator_type::iterator_category> const_iterator;
316        typedef const_iterator iterator;
317#else
318        class const_iterator;
319        typedef const_iterator iterator;
320#endif
321
322        // Element lookup
323        BOOST_UBLAS_INLINE
324        const_iterator find (size_type i) const {
325#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
326            const_subiterator_type it (e_.find (i));
327            return const_iterator (*this, it.index ());
328#else
329            return const_iterator (*this, e_.find (i));
330#endif
331        }
332
333        // Iterator enhances the iterator of the referenced expression
334        // with the unary functor.
335
336#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
337        class const_iterator:
338            public container_const_reference<vector_unary>,
339            public iterator_base_traits<typename E::const_iterator::iterator_category>::template
340                        iterator_base<const_iterator, value_type>::type {
341        public:
342            typedef typename E::const_iterator::iterator_category iterator_category;
343            typedef typename vector_unary::difference_type difference_type;
344            typedef typename vector_unary::value_type value_type;
345            typedef typename vector_unary::const_reference reference;
346            typedef typename vector_unary::const_pointer pointer;
347
348            // Construction and destruction
349            BOOST_UBLAS_INLINE
350            const_iterator ():
351                container_const_reference<self_type> (), it_ () {}
352            BOOST_UBLAS_INLINE
353            const_iterator (const self_type &vu, const const_subiterator_type &it):
354                container_const_reference<self_type> (vu), it_ (it) {}
355
356            // Arithmetic
357            BOOST_UBLAS_INLINE
358            const_iterator &operator ++ () {
359                ++ it_;
360                return *this;
361            }
362            BOOST_UBLAS_INLINE
363            const_iterator &operator -- () {
364                -- it_;
365                return *this;
366            }
367            BOOST_UBLAS_INLINE
368            const_iterator &operator += (difference_type n) {
369                it_ += n;
370                return *this;
371            }
372            BOOST_UBLAS_INLINE
373            const_iterator &operator -= (difference_type n) {
374                it_ -= n;
375                return *this;
376            }
377            BOOST_UBLAS_INLINE
378            difference_type operator - (const const_iterator &it) const {
379                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
380                return it_ - it.it_;
381            }
382
383            // Dereference
384            BOOST_UBLAS_INLINE
385            const_reference operator * () const {
386                return functor_type::apply (*it_);
387            }
388
389            // Index
390            BOOST_UBLAS_INLINE
391            size_type index () const {
392                return it_.index ();
393            }
394
395            // Assignment
396            BOOST_UBLAS_INLINE
397            const_iterator &operator = (const const_iterator &it) {
398                container_const_reference<self_type>::assign (&it ());
399                it_ = it.it_;
400                return *this;
401            }
402
403            // Comparison
404            BOOST_UBLAS_INLINE
405            bool operator == (const const_iterator &it) const {
406                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
407                return it_ == it.it_;
408            }
409            BOOST_UBLAS_INLINE
410            bool operator < (const const_iterator &it) const {
411                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
412                return it_ < it.it_;
413            }
414
415        private:
416            const_subiterator_type it_;
417        };
418#endif
419
420        BOOST_UBLAS_INLINE
421        const_iterator begin () const {
422            return find (0);
423        }
424        BOOST_UBLAS_INLINE
425        const_iterator end () const {
426            return find (size ());
427        }
428
429        // Reverse iterator
430        typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
431
432        BOOST_UBLAS_INLINE
433        const_reverse_iterator rbegin () const {
434            return const_reverse_iterator (end ());
435        }
436        BOOST_UBLAS_INLINE
437        const_reverse_iterator rend () const {
438            return const_reverse_iterator (begin ());
439        }
440
441    private:
442        expression_closure_type e_;
443    };
444
445    template<class E, class F>
446    struct vector_unary_traits {
447        typedef vector_unary<E, F> expression_type;
448//FIXME
449// #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
450        typedef expression_type result_type;
451// #else
452//         typedef typename E::vector_temporary_type result_type;
453// #endif
454    };
455
456    // (- v) [i] = - v [i]
457    template<class E>
458    BOOST_UBLAS_INLINE
459    typename vector_unary_traits<E, scalar_negate<typename E::value_type> >::result_type
460    operator - (const vector_expression<E> &e) {
461        typedef typename vector_unary_traits<E, scalar_negate<typename E::value_type> >::expression_type expression_type;
462        return expression_type (e ());
463    }
464
465    // (conj v) [i] = conj (v [i])
466    template<class E>
467    BOOST_UBLAS_INLINE
468    typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::result_type
469    conj (const vector_expression<E> &e) {
470        typedef typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
471        return expression_type (e ());
472    }
473
474    // (real v) [i] = real (v [i])
475    template<class E>
476    BOOST_UBLAS_INLINE
477    typename vector_unary_traits<E, scalar_real<typename E::value_type> >::result_type
478    real (const vector_expression<E> &e) {
479        typedef typename vector_unary_traits<E, scalar_real<typename E::value_type> >::expression_type expression_type;
480        return expression_type (e ());
481    }
482
483    // (imag v) [i] = imag (v [i])
484    template<class E>
485    BOOST_UBLAS_INLINE
486    typename vector_unary_traits<E, scalar_imag<typename E::value_type> >::result_type
487    imag (const vector_expression<E> &e) {
488        typedef typename vector_unary_traits<E, scalar_imag<typename E::value_type> >::expression_type expression_type;
489        return expression_type (e ());
490    }
491
492    // (trans v) [i] = v [i]
493    template<class E>
494    BOOST_UBLAS_INLINE
495    typename vector_unary_traits<const E, scalar_identity<typename E::value_type> >::result_type
496    trans (const vector_expression<E> &e) {
497        typedef typename vector_unary_traits<const E, scalar_identity<typename E::value_type> >::expression_type expression_type;
498        return expression_type (e ());
499    }
500    template<class E>
501    BOOST_UBLAS_INLINE
502    typename vector_unary_traits<E, scalar_identity<typename E::value_type> >::result_type
503    trans (vector_expression<E> &e) {
504        typedef typename vector_unary_traits<E, scalar_identity<typename E::value_type> >::expression_type expression_type;
505        return expression_type (e ());
506    }
507
508    // (herm v) [i] = conj (v [i])
509    template<class E>
510    BOOST_UBLAS_INLINE
511    typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::result_type
512    herm (const vector_expression<E> &e) {
513        typedef typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
514        return expression_type (e ());
515    }
516
517    template<class E1, class E2, class F>
518    class vector_binary:
519        public vector_expression<vector_binary<E1, E2, F> > {
520
521        typedef E1 expression1_type;
522        typedef E2 expression2_type;
523        typedef F functor_type;
524        typedef typename E1::const_closure_type expression1_closure_type;
525        typedef typename E2::const_closure_type expression2_closure_type;
526        typedef vector_binary<E1, E2, F> self_type;
527    public:
528#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
529        using vector_expression<vector_binary<E1, E2, F> >::operator ();
530#endif
531        typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
532        typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
533        typedef typename F::result_type value_type;
534        typedef value_type const_reference;
535        typedef const_reference reference;
536        typedef const self_type const_closure_type;
537        typedef const_closure_type closure_type;
538        typedef unknown_storage_tag storage_category;
539
540        // Construction and destruction
541        BOOST_UBLAS_INLINE
542        vector_binary (const expression1_type &e1, const expression2_type &e2):
543            e1_ (e1), e2_ (e2) {}
544
545        // Accessors
546        BOOST_UBLAS_INLINE
547        size_type size () const {
548            return BOOST_UBLAS_SAME (e1_.size (), e2_.size ());
549        }
550
551    private:
552        // Accessors
553        BOOST_UBLAS_INLINE
554        const expression1_closure_type &expression1 () const {
555            return e1_;
556        }
557        BOOST_UBLAS_INLINE
558        const expression2_closure_type &expression2 () const {
559            return e2_;
560        }
561
562    public:
563        // Element access
564        BOOST_UBLAS_INLINE
565        const_reference operator () (size_type i) const {
566            return functor_type::apply (e1_ (i), e2_ (i));
567        }
568
569        BOOST_UBLAS_INLINE
570        const_reference operator [] (size_type i) const {
571            return functor_type::apply (e1_ [i], e2_ [i]);
572        }
573
574        // Closure comparison
575        BOOST_UBLAS_INLINE
576        bool same_closure (const vector_binary &vb) const {
577            return (*this).expression1 ().same_closure (vb.expression1 ()) &&
578                   (*this).expression2 ().same_closure (vb.expression2 ());
579        }
580
581        // Iterator types
582    private:
583        typedef typename E1::const_iterator const_subiterator1_type;
584        typedef typename E2::const_iterator const_subiterator2_type;
585        typedef const value_type *const_pointer;
586
587    public:
588#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
589        typedef typename iterator_restrict_traits<typename const_subiterator1_type::iterator_category,
590                                                  typename const_subiterator2_type::iterator_category>::iterator_category iterator_category;
591        typedef indexed_const_iterator<const_closure_type, iterator_category> const_iterator;
592        typedef const_iterator iterator;
593#else
594        class const_iterator;
595        typedef const_iterator iterator;
596#endif
597
598        // Element lookup
599        BOOST_UBLAS_INLINE
600        const_iterator find (size_type i) const {
601            const_subiterator1_type it1 (e1_.find (i));
602            const_subiterator1_type it1_end (e1_.find (size ()));
603            const_subiterator2_type it2 (e2_.find (i));
604            const_subiterator2_type it2_end (e2_.find (size ()));
605            i = (std::min) (it1 != it1_end ? it1.index () : size (),
606                          it2 != it2_end ? it2.index () : size ());
607#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
608            return const_iterator (*this, i);
609#else
610            return const_iterator (*this, i, it1, it1_end, it2, it2_end);
611#endif
612        }
613
614        // Iterator merges the iterators of the referenced expressions and
615        // enhances them with the binary functor.
616
617#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
618        class const_iterator:
619            public container_const_reference<vector_binary>,
620            public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
621                                                                          typename E2::const_iterator::iterator_category>::iterator_category>::template
622                        iterator_base<const_iterator, value_type>::type {
623        public:
624            typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
625                                                      typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
626            typedef typename vector_binary::difference_type difference_type;
627            typedef typename vector_binary::value_type value_type;
628            typedef typename vector_binary::const_reference reference;
629            typedef typename vector_binary::const_pointer pointer;
630
631            // Construction and destruction
632            BOOST_UBLAS_INLINE
633            const_iterator ():
634                container_const_reference<self_type> (), i_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
635            BOOST_UBLAS_INLINE
636            const_iterator (const self_type &vb, size_type i,
637                            const const_subiterator1_type &it1, const const_subiterator1_type &it1_end,
638                            const const_subiterator2_type &it2, const const_subiterator2_type &it2_end):
639                container_const_reference<self_type> (vb), i_ (i), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
640
641            // Dense specializations
642            BOOST_UBLAS_INLINE
643            void increment (dense_random_access_iterator_tag) {
644                ++ i_, ++ it1_, ++ it2_;
645            }
646            BOOST_UBLAS_INLINE
647            void decrement (dense_random_access_iterator_tag) {
648                -- i_, -- it1_, -- it2_;
649            }
650            BOOST_UBLAS_INLINE
651            value_type dereference (dense_random_access_iterator_tag) const {
652                return functor_type::apply (*it1_, *it2_);
653            }
654
655            // Packed specializations
656            BOOST_UBLAS_INLINE
657            void increment (packed_random_access_iterator_tag) {
658                if (it1_ != it1_end_)
659                    if (it1_.index () <= i_)
660                        ++ it1_;
661                if (it2_ != it2_end_)
662                    if (it2_.index () <= i_)
663                        ++ it2_;
664                ++ i_;
665            }
666            BOOST_UBLAS_INLINE
667            void decrement (packed_random_access_iterator_tag) {
668                if (it1_ != it1_end_)
669                    if (i_ <= it1_.index ())
670                        -- it1_;
671                if (it2_ != it2_end_)
672                    if (i_ <= it2_.index ())
673                        -- it2_;
674                -- i_;
675            }
676            BOOST_UBLAS_INLINE
677            value_type dereference (packed_random_access_iterator_tag) const {
678                value_type t1 = value_type/*zero*/();
679                if (it1_ != it1_end_)
680                    if (it1_.index () == i_)
681                        t1 = *it1_;
682                value_type t2 = value_type/*zero*/();
683                if (it2_ != it2_end_)
684                    if (it2_.index () == i_)
685                        t2 = *it2_;
686                return functor_type::apply (t1, t2);
687            }
688
689            // Sparse specializations
690            BOOST_UBLAS_INLINE
691            void increment (sparse_bidirectional_iterator_tag) {
692                size_type index1 = (*this) ().size ();
693                if (it1_ != it1_end_) {
694                    if  (it1_.index () <= i_)
695                        ++ it1_;
696                    if (it1_ != it1_end_)
697                        index1 = it1_.index ();
698                }
699                size_type index2 = (*this) ().size ();
700                if (it2_ != it2_end_) {
701                    if (it2_.index () <= i_)
702                        ++ it2_;
703                    if (it2_ != it2_end_)
704                        index2 = it2_.index ();
705                }
706                i_ = (std::min) (index1, index2);
707            }
708            BOOST_UBLAS_INLINE
709            void decrement (sparse_bidirectional_iterator_tag) {
710                size_type index1 = (*this) ().size ();
711                if (it1_ != it1_end_) {
712                    if (i_ <= it1_.index ())
713                        -- it1_;
714                    if (it1_ != it1_end_)
715                        index1 = it1_.index ();
716                }
717                size_type index2 = (*this) ().size ();
718                if (it2_ != it2_end_) {
719                    if (i_ <= it2_.index ())
720                        -- it2_;
721                    if (it2_ != it2_end_)
722                        index2 = it2_.index ();
723                }
724                i_ = (std::max) (index1, index2);
725            }
726            BOOST_UBLAS_INLINE
727            value_type dereference (sparse_bidirectional_iterator_tag) const {
728                value_type t1 = value_type/*zero*/();
729                if (it1_ != it1_end_)
730                    if (it1_.index () == i_)
731                        t1 = *it1_;
732                value_type t2 = value_type/*zero*/();
733                if (it2_ != it2_end_)
734                    if (it2_.index () == i_)
735                        t2 = *it2_;
736                return functor_type::apply (t1, t2);
737            }
738
739            // Arithmetic
740            BOOST_UBLAS_INLINE
741            const_iterator &operator ++ () {
742                increment (iterator_category ());
743                return *this;
744            }
745            BOOST_UBLAS_INLINE
746            const_iterator &operator -- () {
747                decrement (iterator_category ());
748                return *this;
749            }
750            BOOST_UBLAS_INLINE
751            const_iterator &operator += (difference_type n) {
752                i_ += n, it1_ += n, it2_ += n;
753                return *this;
754            }
755            BOOST_UBLAS_INLINE
756            const_iterator &operator -= (difference_type n) {
757                i_ -= n, it1_ -= n, it2_ -= n;
758                return *this;
759            }
760            BOOST_UBLAS_INLINE
761            difference_type operator - (const const_iterator &it) const {
762                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
763                return index () - it.index ();
764            }
765
766            // Dereference
767            BOOST_UBLAS_INLINE
768            const_reference operator * () const {
769                return dereference (iterator_category ());
770            }
771
772            // Index
773            BOOST_UBLAS_INLINE
774            size_type index () const {
775                return i_;
776            }
777
778            // Assignment
779            BOOST_UBLAS_INLINE
780            const_iterator &operator = (const const_iterator &it) {
781                container_const_reference<self_type>::assign (&it ());
782                i_ = it.i_;
783                it1_ = it.it1_;
784                it1_end_ = it.it1_end_;
785                it2_ = it.it2_;
786                it2_end_ = it.it2_end_;
787                return *this;
788            }
789
790            // Comparison
791            BOOST_UBLAS_INLINE
792            bool operator == (const const_iterator &it) const {
793                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
794                return index () == it.index ();
795            }
796            BOOST_UBLAS_INLINE
797            bool operator < (const const_iterator &it) const {
798                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
799                return index () < it.index ();
800            }
801
802        private:
803            size_type i_;
804            const_subiterator1_type it1_;
805            const_subiterator1_type it1_end_;
806            const_subiterator2_type it2_;
807            const_subiterator2_type it2_end_;
808        };
809#endif
810
811        BOOST_UBLAS_INLINE
812        const_iterator begin () const {
813            return find (0);
814        }
815        BOOST_UBLAS_INLINE
816        const_iterator end () const {
817            return find (size ());
818        }
819
820        // Reverse iterator
821        typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
822
823        BOOST_UBLAS_INLINE
824        const_reverse_iterator rbegin () const {
825            return const_reverse_iterator (end ());
826        }
827        BOOST_UBLAS_INLINE
828        const_reverse_iterator rend () const {
829            return const_reverse_iterator (begin ());
830        }
831
832    private:
833        expression1_closure_type e1_;
834        expression2_closure_type e2_;
835    };
836
837    template<class E1, class E2, class F>
838    struct vector_binary_traits {
839        typedef vector_binary<E1, E2, F> expression_type;
840#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
841        typedef expression_type result_type;
842#else
843        typedef typename E1::vector_temporary_type result_type;
844#endif
845    };
846
847    // (v1 + v2) [i] = v1 [i] + v2 [i]
848    template<class E1, class E2>
849    BOOST_UBLAS_INLINE
850    typename vector_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
851                                                      typename E2::value_type> >::result_type
852    operator + (const vector_expression<E1> &e1,
853                const vector_expression<E2> &e2) {
854        typedef typename vector_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
855                                                                              typename E2::value_type> >::expression_type expression_type;
856        return expression_type (e1 (), e2 ());
857    }
858
859    // (v1 - v2) [i] = v1 [i] - v2 [i]
860    template<class E1, class E2>
861    BOOST_UBLAS_INLINE
862    typename vector_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
863                                                       typename E2::value_type> >::result_type
864    operator - (const vector_expression<E1> &e1,
865                const vector_expression<E2> &e2) {
866        typedef typename vector_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
867                                                                               typename E2::value_type> >::expression_type expression_type;
868        return expression_type (e1 (), e2 ());
869    }
870
871    // (v1 * v2) [i] = v1 [i] * v2 [i]
872    template<class E1, class E2>
873    BOOST_UBLAS_INLINE
874    typename vector_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
875                                                            typename E2::value_type> >::result_type
876    element_prod (const vector_expression<E1> &e1,
877                  const vector_expression<E2> &e2) {
878        typedef typename vector_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
879                                                                                    typename E2::value_type> >::expression_type expression_type;
880        return expression_type (e1 (), e2 ());
881    }
882
883    // (v1 / v2) [i] = v1 [i] / v2 [i]
884    template<class E1, class E2>
885    BOOST_UBLAS_INLINE
886    typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
887                                                         typename E2::value_type> >::result_type
888    element_div (const vector_expression<E1> &e1,
889                 const vector_expression<E2> &e2) {
890        typedef typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
891                                                                                 typename E2::value_type> >::expression_type expression_type;
892        return expression_type (e1 (), e2 ());
893    }
894
895
896    template<class E1, class E2, class F>
897    class vector_binary_scalar1:
898        public vector_expression<vector_binary_scalar1<E1, E2, F> > {
899
900        typedef F functor_type;
901        typedef E1 expression1_type;
902        typedef E2 expression2_type;
903    public:
904        typedef const E1& expression1_closure_type;
905        typedef typename E2::const_closure_type expression2_closure_type;
906    private:
907        typedef vector_binary_scalar1<E1, E2, F> self_type;
908    public:
909#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
910        using vector_expression<vector_binary_scalar1<E1, E2, F> >::operator ();
911#endif
912        typedef typename E2::size_type size_type;
913        typedef typename E2::difference_type difference_type;
914        typedef typename F::result_type value_type;
915        typedef value_type const_reference;
916        typedef const_reference reference;
917        typedef const self_type const_closure_type;
918        typedef const_closure_type closure_type;
919        typedef unknown_storage_tag storage_category;
920
921        // Construction and destruction
922        BOOST_UBLAS_INLINE
923        vector_binary_scalar1 (const expression1_type &e1, const expression2_type &e2):
924            e1_ (e1), e2_ (e2) {}
925
926        // Accessors
927        BOOST_UBLAS_INLINE
928        size_type size () const {
929            return e2_.size ();
930        }
931
932    public:
933        // Element access
934        BOOST_UBLAS_INLINE
935        const_reference operator () (size_type i) const {
936            return functor_type::apply (e1_, e2_ (i));
937        }
938
939        BOOST_UBLAS_INLINE
940        const_reference operator [] (size_type i) const {
941            return functor_type::apply (e1_, e2_ [i]);
942        }
943
944        // Closure comparison
945        BOOST_UBLAS_INLINE
946        bool same_closure (const vector_binary_scalar1 &vbs1) const {
947            return &e1_ == &(vbs1.e1_) &&
948                   (*this).e2_.same_closure (vbs1.e2_);
949        }
950
951        // Iterator types
952    private:
953        typedef expression1_type const_subiterator1_type;
954        typedef typename expression2_type::const_iterator const_subiterator2_type;
955        typedef const value_type *const_pointer;
956
957    public:
958#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
959        typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
960        typedef const_iterator iterator;
961#else
962        class const_iterator;
963        typedef const_iterator iterator;
964#endif
965
966        // Element lookup
967        BOOST_UBLAS_INLINE
968        const_iterator find (size_type i) const {
969#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
970            const_subiterator2_type it (e2_.find (i));
971            return const_iterator (*this, it.index ());
972#else
973            return const_iterator (*this, const_subiterator1_type (e1_), e2_.find (i));
974#endif
975        }
976
977        // Iterator enhances the iterator of the referenced vector expression
978        // with the binary functor.
979
980#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
981        class const_iterator:
982            public container_const_reference<vector_binary_scalar1>,
983            public iterator_base_traits<typename E2::const_iterator::iterator_category>::template
984                        iterator_base<const_iterator, value_type>::type {
985        public:
986            typedef typename E2::const_iterator::iterator_category iterator_category;
987            typedef typename vector_binary_scalar1::difference_type difference_type;
988            typedef typename vector_binary_scalar1::value_type value_type;
989            typedef typename vector_binary_scalar1::const_reference reference;
990            typedef typename vector_binary_scalar1::const_pointer pointer;
991
992            // Construction and destruction
993            BOOST_UBLAS_INLINE
994            const_iterator ():
995                container_const_reference<self_type> (), it1_ (), it2_ () {}
996            BOOST_UBLAS_INLINE
997            const_iterator (const self_type &vbs, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
998                container_const_reference<self_type> (vbs), it1_ (it1), it2_ (it2) {}
999
1000            // Arithmetic
1001            BOOST_UBLAS_INLINE
1002            const_iterator &operator ++ () {
1003                ++ it2_;
1004                return *this;
1005            }
1006            BOOST_UBLAS_INLINE
1007            const_iterator &operator -- () {
1008                -- it2_;
1009                return *this;
1010            }
1011            BOOST_UBLAS_INLINE
1012            const_iterator &operator += (difference_type n) {
1013                it2_ += n;
1014                return *this;
1015            }
1016            BOOST_UBLAS_INLINE
1017            const_iterator &operator -= (difference_type n) {
1018                it2_ -= n;
1019                return *this;
1020            }
1021            BOOST_UBLAS_INLINE
1022            difference_type operator - (const const_iterator &it) const {
1023                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1024                // FIXME we shouldn't compare floats
1025                // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1026                return it2_ - it.it2_;
1027            }
1028
1029            // Dereference
1030            BOOST_UBLAS_INLINE
1031            const_reference operator * () const {
1032                return functor_type::apply (it1_, *it2_);
1033            }
1034
1035            // Index
1036            BOOST_UBLAS_INLINE
1037            size_type index () const {
1038                return it2_.index ();
1039            }
1040
1041            // Assignment
1042            BOOST_UBLAS_INLINE
1043            const_iterator &operator = (const const_iterator &it) {
1044                container_const_reference<self_type>::assign (&it ());
1045                it1_ = it.it1_;
1046                it2_ = it.it2_;
1047                return *this;
1048            }
1049
1050            // Comparison
1051            BOOST_UBLAS_INLINE
1052            bool operator == (const const_iterator &it) const {
1053                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1054                // FIXME we shouldn't compare floats
1055                // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1056                return it2_ == it.it2_;
1057            }
1058            BOOST_UBLAS_INLINE
1059            bool operator < (const const_iterator &it) const {
1060                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1061                // FIXME we shouldn't compare floats
1062                // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1063                return it2_ < it.it2_;
1064            }
1065
1066        private:
1067            const_subiterator1_type it1_;
1068            const_subiterator2_type it2_;
1069        };
1070#endif
1071
1072        BOOST_UBLAS_INLINE
1073        const_iterator begin () const {
1074            return find (0);
1075        }
1076        BOOST_UBLAS_INLINE
1077        const_iterator end () const {
1078            return find (size ());
1079        }
1080
1081        // Reverse iterator
1082        typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1083
1084        BOOST_UBLAS_INLINE
1085        const_reverse_iterator rbegin () const {
1086            return const_reverse_iterator (end ());
1087        }
1088        BOOST_UBLAS_INLINE
1089        const_reverse_iterator rend () const {
1090            return const_reverse_iterator (begin ());
1091        }
1092
1093    private:
1094        expression1_closure_type e1_;
1095        expression2_closure_type e2_;
1096    };
1097
1098    template<class E1, class E2, class F>
1099    struct vector_binary_scalar1_traits {
1100        typedef vector_binary_scalar1<E1, E2, F> expression_type;   // allow E1 to be builtin type
1101#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1102        typedef expression_type result_type;
1103#else
1104        typedef typename E2::vector_temporary_type result_type;
1105#endif
1106    };
1107
1108    // (t * v) [i] = t * v [i]
1109    template<class T1, class E2>
1110    BOOST_UBLAS_INLINE
1111    typename vector_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
1112    operator * (const T1 &e1,
1113                const vector_expression<E2> &e2) {
1114        typedef typename vector_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::expression_type expression_type;
1115        return expression_type (e1, e2 ());
1116    }
1117
1118
1119    template<class E1, class E2, class F>
1120    class vector_binary_scalar2:
1121        public vector_expression<vector_binary_scalar2<E1, E2, F> > {
1122
1123        typedef F functor_type;
1124        typedef E1 expression1_type;
1125        typedef E2 expression2_type;
1126        typedef typename E1::const_closure_type expression1_closure_type;
1127        typedef const E2& expression2_closure_type;
1128        typedef vector_binary_scalar2<E1, E2, F> self_type;
1129    public:
1130#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1131        using vector_expression<vector_binary_scalar2<E1, E2, F> >::operator ();
1132#endif
1133        typedef typename E1::size_type size_type;
1134        typedef typename E1::difference_type difference_type;
1135        typedef typename F::result_type value_type;
1136        typedef value_type const_reference;
1137        typedef const_reference reference;
1138        typedef const self_type const_closure_type;
1139        typedef const_closure_type closure_type;
1140        typedef unknown_storage_tag storage_category;
1141
1142        // Construction and destruction
1143        BOOST_UBLAS_INLINE
1144        vector_binary_scalar2 (const expression1_type &e1, const expression2_type &e2):
1145            e1_ (e1), e2_ (e2) {}
1146
1147        // Accessors
1148        BOOST_UBLAS_INLINE
1149        size_type size () const {
1150            return e1_.size ();
1151        }
1152
1153    public:
1154        // Element access
1155        BOOST_UBLAS_INLINE
1156        const_reference operator () (size_type i) const {
1157            return functor_type::apply (e1_ (i), e2_);
1158        }
1159
1160        BOOST_UBLAS_INLINE
1161        const_reference operator [] (size_type i) const {
1162            return functor_type::apply (e1_ [i], e2_);
1163        }
1164
1165        // Closure comparison
1166        BOOST_UBLAS_INLINE
1167        bool same_closure (const vector_binary_scalar2 &vbs2) const {
1168            return (*this).e1_.same_closure (vbs2.e1_) &&
1169                   &e2_ == &(vbs2.e2_);
1170        }
1171
1172        // Iterator types
1173    private:
1174        typedef typename expression1_type::const_iterator const_subiterator1_type;
1175        typedef expression2_type const_subiterator2_type;
1176        typedef const value_type *const_pointer;
1177
1178    public:
1179#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1180        typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
1181        typedef const_iterator iterator;
1182#else
1183        class const_iterator;
1184        typedef const_iterator iterator;
1185#endif
1186
1187        // Element lookup
1188        BOOST_UBLAS_INLINE
1189        const_iterator find (size_type i) const {
1190#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1191            const_subiterator1_type it (e1_.find (i));
1192            return const_iterator (*this, it.index ());
1193#else
1194            return const_iterator (*this, e1_.find (i), const_subiterator2_type (e2_));
1195#endif
1196        }
1197
1198        // Iterator enhances the iterator of the referenced vector expression
1199        // with the binary functor.
1200
1201#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1202        class const_iterator:
1203            public container_const_reference<vector_binary_scalar2>,
1204            public iterator_base_traits<typename E1::const_iterator::iterator_category>::template
1205                        iterator_base<const_iterator, value_type>::type {
1206        public:
1207            typedef typename E1::const_iterator::iterator_category iterator_category;
1208            typedef typename vector_binary_scalar2::difference_type difference_type;
1209            typedef typename vector_binary_scalar2::value_type value_type;
1210            typedef typename vector_binary_scalar2::const_reference reference;
1211            typedef typename vector_binary_scalar2::const_pointer pointer;
1212
1213            // Construction and destruction
1214            BOOST_UBLAS_INLINE
1215            const_iterator ():
1216                container_const_reference<self_type> (), it1_ (), it2_ () {}
1217            BOOST_UBLAS_INLINE
1218            const_iterator (const self_type &vbs, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
1219                container_const_reference<self_type> (vbs), it1_ (it1), it2_ (it2) {}
1220
1221            // Arithmetic
1222            BOOST_UBLAS_INLINE
1223            const_iterator &operator ++ () {
1224                ++ it1_;
1225                return *this;
1226            }
1227            BOOST_UBLAS_INLINE
1228            const_iterator &operator -- () {
1229                -- it1_;
1230                return *this;
1231            }
1232            BOOST_UBLAS_INLINE
1233            const_iterator &operator += (difference_type n) {
1234                it1_ += n;
1235                return *this;
1236            }
1237            BOOST_UBLAS_INLINE
1238            const_iterator &operator -= (difference_type n) {
1239                it1_ -= n;
1240                return *this;
1241            }
1242            BOOST_UBLAS_INLINE
1243            difference_type operator - (const const_iterator &it) const {
1244                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1245                // FIXME we shouldn't compare floats
1246                // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1247                return it1_ - it.it1_;
1248            }
1249
1250            // Dereference
1251            BOOST_UBLAS_INLINE
1252            const_reference operator * () const {
1253                return functor_type::apply (*it1_, it2_);
1254            }
1255
1256            // Index
1257            BOOST_UBLAS_INLINE
1258            size_type index () const {
1259                return it1_.index ();
1260            }
1261
1262            // Assignment
1263            BOOST_UBLAS_INLINE
1264            const_iterator &operator = (const const_iterator &it) {
1265                container_const_reference<self_type>::assign (&it ());
1266                it1_ = it.it1_;
1267                it2_ = it.it2_;
1268                return *this;
1269            }
1270
1271            // Comparison
1272            BOOST_UBLAS_INLINE
1273            bool operator == (const const_iterator &it) const {
1274                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1275                // FIXME we shouldn't compare floats
1276                // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1277                return it1_ == it.it1_;
1278            }
1279            BOOST_UBLAS_INLINE
1280            bool operator < (const const_iterator &it) const {
1281                BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1282                // FIXME we shouldn't compare floats
1283                // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1284                return it1_ < it.it1_;
1285            }
1286
1287        private:
1288            const_subiterator1_type it1_;
1289            const_subiterator2_type it2_;
1290        };
1291#endif
1292
1293        BOOST_UBLAS_INLINE
1294        const_iterator begin () const {
1295            return find (0);
1296        }
1297        BOOST_UBLAS_INLINE
1298        const_iterator end () const {
1299            return find (size ());
1300        }
1301
1302        // Reverse iterator
1303        typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1304
1305        BOOST_UBLAS_INLINE
1306        const_reverse_iterator rbegin () const {
1307            return const_reverse_iterator (end ());
1308        }
1309        BOOST_UBLAS_INLINE
1310        const_reverse_iterator rend () const {
1311            return const_reverse_iterator (begin ());
1312        }
1313
1314    private:
1315        expression1_closure_type e1_;
1316        expression2_closure_type e2_;
1317    };
1318
1319    template<class E1, class E2, class F>
1320    struct vector_binary_scalar2_traits {
1321        typedef vector_binary_scalar2<E1, E2, F> expression_type;   // allow E2 to be builtin type
1322#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1323        typedef expression_type result_type;
1324#else
1325        typedef typename E1::vector_temporary_type result_type;
1326#endif
1327    };
1328
1329    // (v * t) [i] = v [i] * t
1330    template<class E1, class T2>
1331    BOOST_UBLAS_INLINE
1332    typename vector_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::result_type
1333    operator * (const vector_expression<E1> &e1,
1334                const T2 &e2) {
1335        typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::expression_type expression_type;
1336        return expression_type (e1 (), e2);
1337    }
1338
1339    // (v / t) [i] = v [i] / t
1340    template<class E1, class T2>
1341    BOOST_UBLAS_INLINE
1342    typename vector_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::result_type
1343    operator / (const vector_expression<E1> &e1,
1344                const T2 &e2) {
1345        typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::expression_type expression_type;
1346        return expression_type (e1 (), e2);
1347    }
1348
1349
1350    template<class E, class F>
1351    class vector_scalar_unary:
1352        public scalar_expression<vector_scalar_unary<E, F> > {
1353
1354        typedef E expression_type;
1355        typedef F functor_type;
1356        typedef typename E::const_closure_type expression_closure_type;
1357        typedef typename E::const_iterator::iterator_category iterator_category;
1358        typedef vector_scalar_unary<E, F> self_type;
1359    public:
1360        typedef typename F::size_type size_type;
1361        typedef typename F::difference_type difference_type;
1362        typedef typename F::result_type value_type;
1363        typedef const self_type const_closure_type;
1364        typedef const_closure_type closure_type;
1365        typedef unknown_storage_tag storage_category;
1366
1367        // Construction and destruction
1368        BOOST_UBLAS_INLINE
1369        explicit vector_scalar_unary (const expression_type &e):
1370            e_ (e) {}
1371
1372    private:
1373        // Expression accessors
1374        BOOST_UBLAS_INLINE
1375        const expression_closure_type &expression () const {
1376            return e_;
1377        }
1378
1379    public:
1380        BOOST_UBLAS_INLINE
1381        operator value_type () const {
1382            return evaluate (iterator_category ());
1383        }
1384
1385    private:
1386        // Dense random access specialization
1387        BOOST_UBLAS_INLINE
1388        value_type evaluate (dense_random_access_iterator_tag) const {
1389#ifdef BOOST_UBLAS_USE_INDEXING
1390            return functor_type::apply (e_);
1391#elif BOOST_UBLAS_USE_ITERATING
1392            difference_type size = e_.size ();
1393            return functor_type::apply (size, e_.begin ());
1394#else
1395            difference_type size = e_.size ();
1396            if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
1397                return functor_type::apply (size, e_.begin ());
1398            else
1399                return functor_type::apply (e_);
1400#endif
1401        }
1402
1403        // Packed bidirectional specialization
1404        BOOST_UBLAS_INLINE
1405        value_type evaluate (packed_random_access_iterator_tag) const {
1406            return functor_type::apply (e_.begin (), e_.end ());
1407        }
1408
1409        // Sparse bidirectional specialization
1410        BOOST_UBLAS_INLINE
1411        value_type evaluate (sparse_bidirectional_iterator_tag) const {
1412            return functor_type::apply (e_.begin (), e_.end ());
1413        }
1414
1415    private:
1416        expression_closure_type e_;
1417    };
1418
1419    template<class E, class F>
1420    struct vector_scalar_unary_traits {
1421        typedef vector_scalar_unary<E, F> expression_type;
1422#if !defined (BOOST_UBLAS_SIMPLE_ET_DEBUG) && defined (BOOST_UBLAS_USE_SCALAR_ET)
1423// FIXME don't define USE_SCALAR_ET other then for testing
1424// They do not work for complex types
1425         typedef expression_type result_type;
1426#else
1427         typedef typename F::result_type result_type;
1428#endif
1429    };
1430
1431    // sum v = sum (v [i])
1432    template<class E>
1433    BOOST_UBLAS_INLINE
1434    typename vector_scalar_unary_traits<E, vector_sum<typename E::value_type> >::result_type
1435    sum (const vector_expression<E> &e) {
1436        typedef typename vector_scalar_unary_traits<E, vector_sum<typename E::value_type> >::expression_type expression_type;
1437        return expression_type (e ());
1438    }
1439
1440    // real: norm_1 v = sum (abs (v [i]))
1441    // complex: norm_1 v = sum (abs (real (v [i])) + abs (imag (v [i])))
1442    template<class E>
1443    BOOST_UBLAS_INLINE
1444    typename vector_scalar_unary_traits<E, vector_norm_1<typename E::value_type> >::result_type
1445    norm_1 (const vector_expression<E> &e) {
1446        typedef typename vector_scalar_unary_traits<E, vector_norm_1<typename E::value_type> >::expression_type expression_type;
1447        return expression_type (e ());
1448    }
1449
1450    // real: norm_2 v = sqrt (sum (v [i] * v [i]))
1451    // complex: norm_2 v = sqrt (sum (v [i] * conj (v [i])))
1452    template<class E>
1453    BOOST_UBLAS_INLINE
1454    typename vector_scalar_unary_traits<E, vector_norm_2<typename E::value_type> >::result_type
1455    norm_2 (const vector_expression<E> &e) {
1456        typedef typename vector_scalar_unary_traits<E, vector_norm_2<typename E::value_type> >::expression_type expression_type;
1457        return expression_type (e ());
1458    }
1459
1460    // real: norm_inf v = maximum (abs (v [i]))
1461    // complex: norm_inf v = maximum (maximum (abs (real (v [i])), abs (imag (v [i]))))
1462    template<class E>
1463    BOOST_UBLAS_INLINE
1464    typename vector_scalar_unary_traits<E, vector_norm_inf<typename E::value_type> >::result_type
1465    norm_inf (const vector_expression<E> &e) {
1466        typedef typename vector_scalar_unary_traits<E, vector_norm_inf<typename E::value_type> >::expression_type expression_type;
1467        return expression_type (e ());
1468    }
1469
1470    // real: index_norm_inf v = minimum (i: abs (v [i]) == maximum (abs (v [i])))
1471    template<class E>
1472    BOOST_UBLAS_INLINE
1473    typename vector_scalar_unary_traits<E, vector_index_norm_inf<typename E::value_type> >::result_type
1474    index_norm_inf (const vector_expression<E> &e) {
1475        typedef typename vector_scalar_unary_traits<E, vector_index_norm_inf<typename E::value_type> >::expression_type expression_type;
1476        return expression_type (e ());
1477    }
1478
1479    template<class E1, class E2, class F>
1480    class vector_scalar_binary:
1481        public scalar_expression<vector_scalar_binary<E1, E2, F> > {
1482
1483        typedef E1 expression1_type;
1484        typedef E2 expression2_type;
1485        typedef F functor_type;
1486        typedef typename E1::const_closure_type expression1_closure_type;
1487        typedef typename E2::const_closure_type expression2_closure_type;
1488        typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
1489                                                  typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
1490        typedef vector_scalar_binary<E1, E2, F> self_type;
1491    public:
1492        static const unsigned complexity = 1;
1493        typedef typename F::size_type size_type;
1494        typedef typename F::difference_type difference_type;
1495        typedef typename F::result_type value_type;
1496        typedef const self_type const_closure_type;
1497        typedef const_closure_type closure_type;
1498        typedef unknown_storage_tag storage_category;
1499
1500        // Construction and destruction
1501        BOOST_UBLAS_INLINE
1502        vector_scalar_binary (const expression1_type &e1, const expression2_type  &e2):
1503            e1_ (e1), e2_ (e2) {}
1504
1505    private:
1506        // Accessors
1507        BOOST_UBLAS_INLINE
1508        const expression1_closure_type &expression1 () const {
1509            return e1_;
1510        }
1511        BOOST_UBLAS_INLINE
1512        const expression2_closure_type &expression2 () const {
1513            return e2_;
1514        }
1515
1516    public:
1517        BOOST_UBLAS_INLINE
1518        operator value_type () const {
1519            return evaluate (iterator_category ());
1520        }
1521
1522    private:
1523        // Dense random access specialization
1524        BOOST_UBLAS_INLINE
1525        value_type evaluate (dense_random_access_iterator_tag) const {
1526#ifdef BOOST_UBLAS_USE_INDEXING
1527            return functor_type::apply (e1_, e2_);
1528#elif BOOST_UBLAS_USE_ITERATING
1529            difference_type size = BOOST_UBLAS_SAME (e1_.size (), e2_.size ());
1530            return functor_type::apply (size, e1_.begin (), e2_.begin ());
1531#else
1532            difference_type size = BOOST_UBLAS_SAME (e1_.size (), e2_.size ());
1533            if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
1534                return functor_type::apply (size, e1_.begin (), e2_.begin ());
1535            else
1536                return functor_type::apply (e1_, e2_);
1537#endif
1538        }
1539
1540        // Packed bidirectional specialization
1541        BOOST_UBLAS_INLINE
1542        value_type evaluate (packed_random_access_iterator_tag) const {
1543            return functor_type::apply (e1_.begin (), e1_.end (), e2_.begin (), e2_.end ());
1544        }
1545
1546        // Sparse bidirectional specialization
1547        BOOST_UBLAS_INLINE
1548        value_type evaluate (sparse_bidirectional_iterator_tag) const {
1549            return functor_type::apply (e1_.begin (), e1_.end (), e2_.begin (), e2_.end (), sparse_bidirectional_iterator_tag ());
1550        }
1551
1552    private:
1553        expression1_closure_type e1_;
1554        expression2_closure_type e2_;
1555    };
1556
1557    template<class E1, class E2, class F>
1558    struct vector_scalar_binary_traits {
1559        typedef vector_scalar_binary<E1, E2, F> expression_type;
1560#if !defined (BOOST_UBLAS_SIMPLE_ET_DEBUG) && defined (BOOST_UBLAS_USE_SCALAR_ET)
1561// FIXME don't define USE_SCALAR_ET other then for testing
1562// They do not work for complex types
1563        typedef expression_type result_type;
1564#else
1565        typedef typename F::result_type result_type;
1566#endif
1567    };
1568
1569    // inner_prod (v1, v2) = sum (v1 [i] * v2 [i])
1570    template<class E1, class E2>
1571    BOOST_UBLAS_INLINE
1572    typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<typename E1::value_type,
1573                                                                   typename E2::value_type,
1574                                                                   typename promote_traits<typename E1::value_type,
1575                                                                                           typename E2::value_type>::promote_type> >::result_type
1576    inner_prod (const vector_expression<E1> &e1,
1577                const vector_expression<E2> &e2) {
1578        typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<typename E1::value_type,
1579                                                                                           typename E2::value_type,
1580                                                                                           typename promote_traits<typename E1::value_type,
1581                                                                                                                               typename E2::value_type>::promote_type> >::expression_type expression_type;
1582        return expression_type (e1 (), e2 ());
1583    }
1584
1585    template<class E1, class E2>
1586    BOOST_UBLAS_INLINE
1587    typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<typename E1::value_type,
1588                                                                   typename E2::value_type,
1589                                                                   typename type_traits<typename promote_traits<typename E1::value_type,
1590                                                                                                                typename E2::value_type>::promote_type>::precision_type> >::result_type
1591    prec_inner_prod (const vector_expression<E1> &e1,
1592                     const vector_expression<E2> &e2) {
1593        typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<typename E1::value_type,
1594                                                                                           typename E2::value_type,
1595                                                                                           typename type_traits<typename promote_traits<typename E1::value_type,
1596                                                                                                                                                                typename E2::value_type>::promote_type>::precision_type> >::expression_type expression_type;
1597        return expression_type (e1 (), e2 ());
1598    }
1599
1600}}}
1601
1602#endif
Note: See TracBrowser for help on using the repository browser.