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