Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

matrix_expression.hpp @ 857

Revision 857, 187.3 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_MATRIX_EXPRESSION_
18	#define _BOOST_UBLAS_MATRIX_EXPRESSION_
19
20	#include <boost/numeric/ublas/vector_expression.hpp>
21
22	// Expression templates based on ideas of Todd Veldhuizen and Geoffrey Furnish
23	// Iterators based on ideas of Jeremy Siek
24	//
25	// Classes that model the Matrix Expression concept
26
27	namespace boost { namespace numeric { namespace ublas {
28
29	template<class E>
30	class matrix_reference:
31	public matrix_expression<matrix_reference<E> > {
32
33	typedef matrix_reference<E> self_type;
34	public:
35	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
36	using matrix_expression<self_type>::operator ();
37	#endif
38	typedef typename E::size_type size_type;
39	typedef typename E::difference_type difference_type;
40	typedef typename E::value_type value_type;
41	typedef typename E::const_reference const_reference;
42	typedef typename boost::mpl::if_<boost::is_const<E>,
43	typename E::const_reference,
44	typename E::reference>::type reference;
45	typedef E refered_type;
46	typedef const self_type const_closure_type;
47	typedef const_closure_type closure_type;
48	typedef typename E::orientation_category orientation_category;
49	typedef typename E::storage_category storage_category;
50
51	// Construction and destruction
52	BOOST_UBLAS_INLINE
53	explicit matrix_reference (refered_type &e):
54	e_ (e) {}
55
56	// Accessors
57	BOOST_UBLAS_INLINE
58	size_type size1 () const {
59	return e_.size1 ();
60	}
61	BOOST_UBLAS_INLINE
62	size_type size2 () const {
63	return e_.size2 ();
64	}
65
66	public:
67	// Expression accessors - const correct
68	BOOST_UBLAS_INLINE
69	const refered_type &expression () const {
70	return e_;
71	}
72	BOOST_UBLAS_INLINE
73	refered_type &expression () {
74	return e_;
75	}
76
77	public:
78	// Element access
79	#ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
80	BOOST_UBLAS_INLINE
81	const_reference operator () (size_type i, size_type j) const {
82	return expression () (i, j);
83	}
84	BOOST_UBLAS_INLINE
85	reference operator () (size_type i, size_type j) {
86	return expression () (i, j);
87	}
88	#else
89	BOOST_UBLAS_INLINE
90	reference operator () (size_type i, size_type j) const {
91	return expression () (i, j);
92	}
93	#endif
94
95	// Assignment
96	BOOST_UBLAS_INLINE
97	matrix_reference &operator = (const matrix_reference &m) {
98	expression ().operator = (m);
99	return *this;
100	}
101	template<class AE>
102	BOOST_UBLAS_INLINE
103	matrix_reference &operator = (const matrix_expression<AE> &ae) {
104	expression ().operator = (ae);
105	return *this;
106	}
107	template<class AE>
108	BOOST_UBLAS_INLINE
109	matrix_reference &assign (const matrix_expression<AE> &ae) {
110	expression ().assign (ae);
111	return *this;
112	}
113	template<class AE>
114	BOOST_UBLAS_INLINE
115	matrix_reference &operator += (const matrix_expression<AE> &ae) {
116	expression ().operator += (ae);
117	return *this;
118	}
119	template<class AE>
120	BOOST_UBLAS_INLINE
121	matrix_reference &plus_assign (const matrix_expression<AE> &ae) {
122	expression ().plus_assign (ae);
123	return *this;
124	}
125	template<class AE>
126	BOOST_UBLAS_INLINE
127	matrix_reference &operator -= (const matrix_expression<AE> &ae) {
128	expression ().operator -= (ae);
129	return *this;
130	}
131	template<class AE>
132	BOOST_UBLAS_INLINE
133	matrix_reference &minus_assign (const matrix_expression<AE> &ae) {
134	expression ().minus_assign (ae);
135	return *this;
136	}
137	template<class AT>
138	BOOST_UBLAS_INLINE
139	matrix_reference &operator *= (const AT &at) {
140	expression ().operator *= (at);
141	return *this;
142	}
143	template<class AT>
144	BOOST_UBLAS_INLINE
145	matrix_reference &operator /= (const AT &at) {
146	expression ().operator /= (at);
147	return *this;
148	}
149
150	// Swapping
151	BOOST_UBLAS_INLINE
152	void swap (matrix_reference &m) {
153	expression ().swap (m.expression ());
154	}
155
156	// Closure comparison
157	BOOST_UBLAS_INLINE
158	bool same_closure (const matrix_reference &mr) const {
159	return &(*this).e_ == &mr.e_;
160	}
161
162	// Iterator types
163	typedef typename E::const_iterator1 const_iterator1;
164	typedef typename boost::mpl::if_<boost::is_const<E>,
165	typename E::const_iterator1,
166	typename E::iterator1>::type iterator1;
167	typedef typename E::const_iterator2 const_iterator2;
168	typedef typename boost::mpl::if_<boost::is_const<E>,
169	typename E::const_iterator2,
170	typename E::iterator2>::type iterator2;
171
172	// Element lookup
173	BOOST_UBLAS_INLINE
174	const_iterator1 find1 (int rank, size_type i, size_type j) const {
175	return expression ().find1 (rank, i, j);
176	}
177	BOOST_UBLAS_INLINE
178	iterator1 find1 (int rank, size_type i, size_type j) {
179	return expression ().find1 (rank, i, j);
180	}
181	BOOST_UBLAS_INLINE
182	const_iterator2 find2 (int rank, size_type i, size_type j) const {
183	return expression ().find2 (rank, i, j);
184	}
185	BOOST_UBLAS_INLINE
186	iterator2 find2 (int rank, size_type i, size_type j) {
187	return expression ().find2 (rank, i, j);
188	}
189
190	// Iterators are the iterators of the referenced expression.
191
192	BOOST_UBLAS_INLINE
193	const_iterator1 begin1 () const {
194	return expression ().begin1 ();
195	}
196	BOOST_UBLAS_INLINE
197	const_iterator1 end1 () const {
198	return expression ().end1 ();
199	}
200
201	BOOST_UBLAS_INLINE
202	iterator1 begin1 () {
203	return expression ().begin1 ();
204	}
205	BOOST_UBLAS_INLINE
206	iterator1 end1 () {
207	return expression ().end1 ();
208	}
209
210	BOOST_UBLAS_INLINE
211	const_iterator2 begin2 () const {
212	return expression ().begin2 ();
213	}
214	BOOST_UBLAS_INLINE
215	const_iterator2 end2 () const {
216	return expression ().end2 ();
217	}
218
219	BOOST_UBLAS_INLINE
220	iterator2 begin2 () {
221	return expression ().begin2 ();
222	}
223	BOOST_UBLAS_INLINE
224	iterator2 end2 () {
225	return expression ().end2 ();
226	}
227
228	// Reverse iterators
229	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
230	typedef reverse_iterator_base1<iterator1> reverse_iterator1;
231
232	BOOST_UBLAS_INLINE
233	const_reverse_iterator1 rbegin1 () const {
234	return const_reverse_iterator1 (end1 ());
235	}
236	BOOST_UBLAS_INLINE
237	const_reverse_iterator1 rend1 () const {
238	return const_reverse_iterator1 (begin1 ());
239	}
240
241	BOOST_UBLAS_INLINE
242	reverse_iterator1 rbegin1 () {
243	return reverse_iterator1 (end1 ());
244	}
245	BOOST_UBLAS_INLINE
246	reverse_iterator1 rend1 () {
247	return reverse_iterator1 (begin1 ());
248	}
249
250	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
251	typedef reverse_iterator_base2<iterator2> reverse_iterator2;
252
253	BOOST_UBLAS_INLINE
254	const_reverse_iterator2 rbegin2 () const {
255	return const_reverse_iterator2 (end2 ());
256	}
257	BOOST_UBLAS_INLINE
258	const_reverse_iterator2 rend2 () const {
259	return const_reverse_iterator2 (begin2 ());
260	}
261
262	BOOST_UBLAS_INLINE
263	reverse_iterator2 rbegin2 () {
264	return reverse_iterator2 (end2 ());
265	}
266	BOOST_UBLAS_INLINE
267	reverse_iterator2 rend2 () {
268	return reverse_iterator2 (begin2 ());
269	}
270
271	private:
272	refered_type &e_;
273	};
274
275
276	template<class E1, class E2, class F>
277	class vector_matrix_binary:
278	public matrix_expression<vector_matrix_binary<E1, E2, F> > {
279
280	typedef E1 expression1_type;
281	typedef E2 expression2_type;
282	public:
283	typedef typename E1::const_closure_type expression1_closure_type;
284	typedef typename E2::const_closure_type expression2_closure_type;
285	private:
286	typedef vector_matrix_binary<E1, E2, F> self_type;
287	public:
288	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
289	using matrix_expression<self_type>::operator ();
290	#endif
291	typedef F functor_type;
292	typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
293	typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
294	typedef typename F::result_type value_type;
295	typedef value_type const_reference;
296	typedef const_reference reference;
297	typedef const self_type const_closure_type;
298	typedef const_closure_type closure_type;
299	typedef unknown_orientation_tag orientation_category;
300	typedef unknown_storage_tag storage_category;
301
302	// Construction and destruction
303	BOOST_UBLAS_INLINE
304	vector_matrix_binary (const expression1_type &e1, const expression2_type &e2):
305	e1_ (e1), e2_ (e2) {}
306
307	// Accessors
308	BOOST_UBLAS_INLINE
309	size_type size1 () const {
310	return e1_.size ();
311	}
312	BOOST_UBLAS_INLINE
313	size_type size2 () const {
314	return e2_.size ();
315	}
316
317	public:
318	// Expression accessors
319	BOOST_UBLAS_INLINE
320	const expression1_closure_type &expression1 () const {
321	return e1_;
322	}
323	BOOST_UBLAS_INLINE
324	const expression2_closure_type &expression2 () const {
325	return e2_;
326	}
327
328	public:
329	// Element access
330	BOOST_UBLAS_INLINE
331	const_reference operator () (size_type i, size_type j) const {
332	return functor_type::apply (e1_ (i), e2_ (j));
333	}
334
335	// Closure comparison
336	BOOST_UBLAS_INLINE
337	bool same_closure (const vector_matrix_binary &vmb) const {
338	return (*this).expression1 ().same_closure (vmb.expression1 ()) &&
339	(*this).expression2 ().same_closure (vmb.expression2 ());
340	}
341
342	// Iterator types
343	private:
344	typedef typename E1::const_iterator const_subiterator1_type;
345	typedef typename E2::const_iterator const_subiterator2_type;
346	typedef const value_type *const_pointer;
347
348	public:
349	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
350	typedef typename iterator_restrict_traits<typename const_subiterator1_type::iterator_category,
351	typename const_subiterator2_type::iterator_category>::iterator_category iterator_category;
352	typedef indexed_const_iterator1<const_closure_type, iterator_category> const_iterator1;
353	typedef const_iterator1 iterator1;
354	typedef indexed_const_iterator2<const_closure_type, iterator_category> const_iterator2;
355	typedef const_iterator2 iterator2;
356	#else
357	class const_iterator1;
358	typedef const_iterator1 iterator1;
359	class const_iterator2;
360	typedef const_iterator2 iterator2;
361	#endif
362	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
363	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
364
365	// Element lookup
366	BOOST_UBLAS_INLINE
367	const_iterator1 find1 (int rank, size_type i, size_type j) const {
368	const_subiterator1_type it1 (e1_.find (i));
369	const_subiterator1_type it1_end (e1_.find (size1 ()));
370	const_subiterator2_type it2 (e2_.find (j));
371	const_subiterator2_type it2_end (e2_.find (size2 ()));
372	if (it2 == it2_end \|\| (rank == 1 && (it2.index () != j \|\| it2 == value_type/zero*/()))) {
373	it1 = it1_end;
374	it2 = it2_end;
375	}
376	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
377	return const_iterator1 (*this, it1.index (), it2.index ());
378	#else
379	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
380	return const_iterator1 (this, it1, it2, it2 != it2_end ? it2 : value_type/zero/());
381	#else
382	return const_iterator1 (*this, it1, it2);
383	#endif
384	#endif
385	}
386	BOOST_UBLAS_INLINE
387	const_iterator2 find2 (int rank, size_type i, size_type j) const {
388	const_subiterator2_type it2 (e2_.find (j));
389	const_subiterator2_type it2_end (e2_.find (size2 ()));
390	const_subiterator1_type it1 (e1_.find (i));
391	const_subiterator1_type it1_end (e1_.find (size1 ()));
392	if (it1 == it1_end \|\| (rank == 1 && (it1.index () != i \|\| it1 == value_type/zero*/()))) {
393	it2 = it2_end;
394	it1 = it1_end;
395	}
396	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
397	return const_iterator2 (*this, it1.index (), it2.index ());
398	#else
399	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
400	return const_iterator2 (this, it1, it2, it1 != it1_end ? it1 : value_type/zero/());
401	#else
402	return const_iterator2 (*this, it1, it2);
403	#endif
404	#endif
405	}
406
407	// Iterators enhance the iterators of the referenced expressions
408	// with the binary functor.
409
410	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
411	class const_iterator1:
412	public container_const_reference<vector_matrix_binary>,
413	public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
414	typename E2::const_iterator::iterator_category>::iterator_category>::template
415	iterator_base<const_iterator1, value_type>::type {
416	public:
417	typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
418	typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
419	typedef typename vector_matrix_binary::difference_type difference_type;
420	typedef typename vector_matrix_binary::value_type value_type;
421	typedef typename vector_matrix_binary::const_reference reference;
422	typedef typename vector_matrix_binary::const_pointer pointer;
423
424	typedef const_iterator2 dual_iterator_type;
425	typedef const_reverse_iterator2 dual_reverse_iterator_type;
426
427	// Construction and destruction
428	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
429	BOOST_UBLAS_INLINE
430	const_iterator1 ():
431	container_const_reference<self_type> (), it1_ (), it2_ (), t2_ () {}
432	BOOST_UBLAS_INLINE
433	const_iterator1 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2, value_type t2):
434	container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2), t2_ (t2) {}
435	#else
436	BOOST_UBLAS_INLINE
437	const_iterator1 ():
438	container_const_reference<self_type> (), it1_ (), it2_ () {}
439	BOOST_UBLAS_INLINE
440	const_iterator1 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
441	container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2) {}
442	#endif
443
444	// Arithmetic
445	BOOST_UBLAS_INLINE
446	const_iterator1 &operator ++ () {
447	++ it1_;
448	return *this;
449	}
450	BOOST_UBLAS_INLINE
451	const_iterator1 &operator -- () {
452	-- it1_;
453	return *this;
454	}
455	BOOST_UBLAS_INLINE
456	const_iterator1 &operator += (difference_type n) {
457	it1_ += n;
458	return *this;
459	}
460	BOOST_UBLAS_INLINE
461	const_iterator1 &operator -= (difference_type n) {
462	it1_ -= n;
463	return *this;
464	}
465	BOOST_UBLAS_INLINE
466	difference_type operator - (const const_iterator1 &it) const {
467	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
468	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
469	return it1_ - it.it1_;
470	}
471
472	// Dereference
473	BOOST_UBLAS_INLINE
474	const_reference operator * () const {
475	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
476	return functor_type::apply (*it1_, t2_);
477	#else
478	return functor_type::apply (it1_, it2_);
479	#endif
480	}
481
482	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
483	BOOST_UBLAS_INLINE
484	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
485	typename self_type::
486	#endif
487	const_iterator2 begin () const {
488	return (*this) ().find2 (1, index1 (), 0);
489	}
490	BOOST_UBLAS_INLINE
491	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
492	typename self_type::
493	#endif
494	const_iterator2 end () const {
495	return (this) ().find2 (1, index1 (), (this) ().size2 ());
496	}
497	BOOST_UBLAS_INLINE
498	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
499	typename self_type::
500	#endif
501	const_reverse_iterator2 rbegin () const {
502	return const_reverse_iterator2 (end ());
503	}
504	BOOST_UBLAS_INLINE
505	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
506	typename self_type::
507	#endif
508	const_reverse_iterator2 rend () const {
509	return const_reverse_iterator2 (begin ());
510	}
511	#endif
512
513	// Indices
514	BOOST_UBLAS_INLINE
515	size_type index1 () const {
516	return it1_.index ();
517	}
518	BOOST_UBLAS_INLINE
519	size_type index2 () const {
520	return it2_.index ();
521	}
522
523	// Assignment
524	BOOST_UBLAS_INLINE
525	const_iterator1 &operator = (const const_iterator1 &it) {
526	container_const_reference<self_type>::assign (&it ());
527	it1_ = it.it1_;
528	it2_ = it.it2_;
529	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
530	t2_ = it.t2_;
531	#endif
532	return *this;
533	}
534
535	// Comparison
536	BOOST_UBLAS_INLINE
537	bool operator == (const const_iterator1 &it) const {
538	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
539	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
540	return it1_ == it.it1_;
541	}
542	BOOST_UBLAS_INLINE
543	bool operator < (const const_iterator1 &it) const {
544	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
545	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
546	return it1_ < it.it1_;
547	}
548
549	private:
550	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
551	const_subiterator1_type it1_;
552	// Mutable due to assignment
553	/* const */ const_subiterator2_type it2_;
554	value_type t2_;
555	#else
556	const_subiterator1_type it1_;
557	const_subiterator2_type it2_;
558	#endif
559	};
560	#endif
561
562	BOOST_UBLAS_INLINE
563	const_iterator1 begin1 () const {
564	return find1 (0, 0, 0);
565	}
566	BOOST_UBLAS_INLINE
567	const_iterator1 end1 () const {
568	return find1 (0, size1 (), 0);
569	}
570
571	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
572	class const_iterator2:
573	public container_const_reference<vector_matrix_binary>,
574	public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
575	typename E2::const_iterator::iterator_category>::iterator_category>::template
576	iterator_base<const_iterator2, value_type>::type {
577	public:
578	typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
579	typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
580	typedef typename vector_matrix_binary::difference_type difference_type;
581	typedef typename vector_matrix_binary::value_type value_type;
582	typedef typename vector_matrix_binary::const_reference reference;
583	typedef typename vector_matrix_binary::const_pointer pointer;
584
585	typedef const_iterator1 dual_iterator_type;
586	typedef const_reverse_iterator1 dual_reverse_iterator_type;
587
588	// Construction and destruction
589	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
590	BOOST_UBLAS_INLINE
591	const_iterator2 ():
592	container_const_reference<self_type> (), it1_ (), it2_ (), t1_ () {}
593	BOOST_UBLAS_INLINE
594	const_iterator2 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2, value_type t1):
595	container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2), t1_ (t1) {}
596	#else
597	BOOST_UBLAS_INLINE
598	const_iterator2 ():
599	container_const_reference<self_type> (), it1_ (), it2_ () {}
600	BOOST_UBLAS_INLINE
601	const_iterator2 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
602	container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2) {}
603	#endif
604
605	// Arithmetic
606	BOOST_UBLAS_INLINE
607	const_iterator2 &operator ++ () {
608	++ it2_;
609	return *this;
610	}
611	BOOST_UBLAS_INLINE
612	const_iterator2 &operator -- () {
613	-- it2_;
614	return *this;
615	}
616	BOOST_UBLAS_INLINE
617	const_iterator2 &operator += (difference_type n) {
618	it2_ += n;
619	return *this;
620	}
621	BOOST_UBLAS_INLINE
622	const_iterator2 &operator -= (difference_type n) {
623	it2_ -= n;
624	return *this;
625	}
626	BOOST_UBLAS_INLINE
627	difference_type operator - (const const_iterator2 &it) const {
628	BOOST_UBLAS_CHECK ((*this) ().same_closure(it ()), external_logic ());
629	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
630	return it2_ - it.it2_;
631	}
632
633	// Dereference
634	BOOST_UBLAS_INLINE
635	const_reference operator * () const {
636	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
637	return functor_type::apply (t1_, *it2_);
638	#else
639	return functor_type::apply (it1_, it2_);
640	#endif
641	}
642
643	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
644	BOOST_UBLAS_INLINE
645	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
646	typename self_type::
647	#endif
648	const_iterator1 begin () const {
649	return (*this) ().find1 (1, 0, index2 ());
650	}
651	BOOST_UBLAS_INLINE
652	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
653	typename self_type::
654	#endif
655	const_iterator1 end () const {
656	return (this) ().find1 (1, (this) ().size1 (), index2 ());
657	}
658	BOOST_UBLAS_INLINE
659	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
660	typename self_type::
661	#endif
662	const_reverse_iterator1 rbegin () const {
663	return const_reverse_iterator1 (end ());
664	}
665	BOOST_UBLAS_INLINE
666	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
667	typename self_type::
668	#endif
669	const_reverse_iterator1 rend () const {
670	return const_reverse_iterator1 (begin ());
671	}
672	#endif
673
674	// Indices
675	BOOST_UBLAS_INLINE
676	size_type index1 () const {
677	return it1_.index ();
678	}
679	BOOST_UBLAS_INLINE
680	size_type index2 () const {
681	return it2_.index ();
682	}
683
684	// Assignment
685	BOOST_UBLAS_INLINE
686	const_iterator2 &operator = (const const_iterator2 &it) {
687	container_const_reference<self_type>::assign (&it ());
688	it1_ = it.it1_;
689	it2_ = it.it2_;
690	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
691	t1_ = it.t1_;
692	#endif
693	return *this;
694	}
695
696	// Comparison
697	BOOST_UBLAS_INLINE
698	bool operator == (const const_iterator2 &it) const {
699	BOOST_UBLAS_CHECK ((*this) ().same_closure( it ()), external_logic ());
700	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
701	return it2_ == it.it2_;
702	}
703	BOOST_UBLAS_INLINE
704	bool operator < (const const_iterator2 &it) const {
705	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
706	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
707	return it2_ < it.it2_;
708	}
709
710	private:
711	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
712	// Mutable due to assignment
713	/* const */ const_subiterator1_type it1_;
714	const_subiterator2_type it2_;
715	value_type t1_;
716	#else
717	const_subiterator1_type it1_;
718	const_subiterator2_type it2_;
719	#endif
720	};
721	#endif
722
723	BOOST_UBLAS_INLINE
724	const_iterator2 begin2 () const {
725	return find2 (0, 0, 0);
726	}
727	BOOST_UBLAS_INLINE
728	const_iterator2 end2 () const {
729	return find2 (0, 0, size2 ());
730	}
731
732	// Reverse iterators
733
734	BOOST_UBLAS_INLINE
735	const_reverse_iterator1 rbegin1 () const {
736	return const_reverse_iterator1 (end1 ());
737	}
738	BOOST_UBLAS_INLINE
739	const_reverse_iterator1 rend1 () const {
740	return const_reverse_iterator1 (begin1 ());
741	}
742
743	BOOST_UBLAS_INLINE
744	const_reverse_iterator2 rbegin2 () const {
745	return const_reverse_iterator2 (end2 ());
746	}
747	BOOST_UBLAS_INLINE
748	const_reverse_iterator2 rend2 () const {
749	return const_reverse_iterator2 (begin2 ());
750	}
751
752	private:
753	expression1_closure_type e1_;
754	expression2_closure_type e2_;
755	};
756
757	template<class E1, class E2, class F>
758	struct vector_matrix_binary_traits {
759	typedef vector_matrix_binary<E1, E2, F> expression_type;
760	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
761	typedef expression_type result_type;
762	#else
763	// ISSUE matrix is arbitary temporary type
764	typedef matrix<typename F::value_type> result_type;
765	#endif
766	};
767
768	// (outer_prod (v1, v2)) [i] [j] = v1 [i] * v2 [j]
769	template<class E1, class E2>
770	BOOST_UBLAS_INLINE
771	typename vector_matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type, typename E2::value_type> >::result_type
772	outer_prod (const vector_expression<E1> &e1,
773	const vector_expression<E2> &e2) {
774	BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
775	typedef typename vector_matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type, typename E2::value_type> >::expression_type expression_type;
776	return expression_type (e1 (), e2 ());
777	}
778
779	template<class E, class F>
780	class matrix_unary1:
781	public matrix_expression<matrix_unary1<E, F> > {
782
783	typedef E expression_type;
784	typedef F functor_type;
785	public:
786	typedef typename E::const_closure_type expression_closure_type;
787	private:
788	typedef matrix_unary1<E, F> self_type;
789	public:
790	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
791	using matrix_expression<self_type>::operator ();
792	#endif
793	typedef typename E::size_type size_type;
794	typedef typename E::difference_type difference_type;
795	typedef typename F::result_type value_type;
796	typedef value_type const_reference;
797	typedef const_reference reference;
798	typedef const self_type const_closure_type;
799	typedef const_closure_type closure_type;
800	typedef typename E::orientation_category orientation_category;
801	typedef unknown_storage_tag storage_category;
802
803	// Construction and destruction
804	BOOST_UBLAS_INLINE
805	explicit matrix_unary1 (const expression_type &e):
806	e_ (e) {}
807
808	// Accessors
809	BOOST_UBLAS_INLINE
810	size_type size1 () const {
811	return e_.size1 ();
812	}
813	BOOST_UBLAS_INLINE
814	size_type size2 () const {
815	return e_.size2 ();
816	}
817
818	public:
819	// Expression accessors
820	BOOST_UBLAS_INLINE
821	const expression_closure_type &expression () const {
822	return e_;
823	}
824
825	public:
826	// Element access
827	BOOST_UBLAS_INLINE
828	const_reference operator () (size_type i, size_type j) const {
829	return functor_type::apply (e_ (i, j));
830	}
831
832	// Closure comparison
833	BOOST_UBLAS_INLINE
834	bool same_closure (const matrix_unary1 &mu1) const {
835	return (*this).expression ().same_closure (mu1.expression ());
836	}
837
838	// Iterator types
839	private:
840	typedef typename E::const_iterator1 const_subiterator1_type;
841	typedef typename E::const_iterator2 const_subiterator2_type;
842	typedef const value_type *const_pointer;
843
844	public:
845	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
846	typedef indexed_const_iterator1<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator1;
847	typedef const_iterator1 iterator1;
848	typedef indexed_const_iterator2<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator2;
849	typedef const_iterator2 iterator2;
850	#else
851	class const_iterator1;
852	typedef const_iterator1 iterator1;
853	class const_iterator2;
854	typedef const_iterator2 iterator2;
855	#endif
856	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
857	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
858
859	// Element lookup
860	BOOST_UBLAS_INLINE
861	const_iterator1 find1 (int rank, size_type i, size_type j) const {
862	const_subiterator1_type it1 (e_.find1 (rank, i, j));
863	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
864	return const_iterator1 (*this, it1.index1 (), it1.index2 ());
865	#else
866	return const_iterator1 (*this, it1);
867	#endif
868	}
869	BOOST_UBLAS_INLINE
870	const_iterator2 find2 (int rank, size_type i, size_type j) const {
871	const_subiterator2_type it2 (e_.find2 (rank, i, j));
872	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
873	return const_iterator2 (*this, it2.index1 (), it2.index2 ());
874	#else
875	return const_iterator2 (*this, it2);
876	#endif
877	}
878
879	// Iterators enhance the iterators of the referenced expression
880	// with the unary functor.
881
882	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
883	class const_iterator1:
884	public container_const_reference<matrix_unary1>,
885	public iterator_base_traits<typename E::const_iterator1::iterator_category>::template
886	iterator_base<const_iterator1, value_type>::type {
887	public:
888	typedef typename E::const_iterator1::iterator_category iterator_category;
889	typedef typename matrix_unary1::difference_type difference_type;
890	typedef typename matrix_unary1::value_type value_type;
891	typedef typename matrix_unary1::const_reference reference;
892	typedef typename matrix_unary1::const_pointer pointer;
893
894	typedef const_iterator2 dual_iterator_type;
895	typedef const_reverse_iterator2 dual_reverse_iterator_type;
896
897	// Construction and destruction
898	BOOST_UBLAS_INLINE
899	const_iterator1 ():
900	container_const_reference<self_type> (), it_ () {}
901	BOOST_UBLAS_INLINE
902	const_iterator1 (const self_type &mu, const const_subiterator1_type &it):
903	container_const_reference<self_type> (mu), it_ (it) {}
904
905	// Arithmetic
906	BOOST_UBLAS_INLINE
907	const_iterator1 &operator ++ () {
908	++ it_;
909	return *this;
910	}
911	BOOST_UBLAS_INLINE
912	const_iterator1 &operator -- () {
913	-- it_;
914	return *this;
915	}
916	BOOST_UBLAS_INLINE
917	const_iterator1 &operator += (difference_type n) {
918	it_ += n;
919	return *this;
920	}
921	BOOST_UBLAS_INLINE
922	const_iterator1 &operator -= (difference_type n) {
923	it_ -= n;
924	return *this;
925	}
926	BOOST_UBLAS_INLINE
927	difference_type operator - (const const_iterator1 &it) const {
928	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
929	return it_ - it.it_;
930	}
931
932	// Dereference
933	BOOST_UBLAS_INLINE
934	const_reference operator * () const {
935	return functor_type::apply (*it_);
936	}
937
938	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
939	BOOST_UBLAS_INLINE
940	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
941	typename self_type::
942	#endif
943	const_iterator2 begin () const {
944	return (*this) ().find2 (1, index1 (), 0);
945	}
946	BOOST_UBLAS_INLINE
947	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
948	typename self_type::
949	#endif
950	const_iterator2 end () const {
951	return (this) ().find2 (1, index1 (), (this) ().size2 ());
952	}
953	BOOST_UBLAS_INLINE
954	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
955	typename self_type::
956	#endif
957	const_reverse_iterator2 rbegin () const {
958	return const_reverse_iterator2 (end ());
959	}
960	BOOST_UBLAS_INLINE
961	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
962	typename self_type::
963	#endif
964	const_reverse_iterator2 rend () const {
965	return const_reverse_iterator2 (begin ());
966	}
967	#endif
968
969	// Indices
970	BOOST_UBLAS_INLINE
971	size_type index1 () const {
972	return it_.index1 ();
973	}
974	BOOST_UBLAS_INLINE
975	size_type index2 () const {
976	return it_.index2 ();
977	}
978
979	// Assignment
980	BOOST_UBLAS_INLINE
981	const_iterator1 &operator = (const const_iterator1 &it) {
982	container_const_reference<self_type>::assign (&it ());
983	it_ = it.it_;
984	return *this;
985	}
986
987	// Comparison
988	BOOST_UBLAS_INLINE
989	bool operator == (const const_iterator1 &it) const {
990	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
991	return it_ == it.it_;
992	}
993	BOOST_UBLAS_INLINE
994	bool operator < (const const_iterator1 &it) const {
995	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
996	return it_ < it.it_;
997	}
998
999	private:
1000	const_subiterator1_type it_;
1001	};
1002	#endif
1003
1004	BOOST_UBLAS_INLINE
1005	const_iterator1 begin1 () const {
1006	return find1 (0, 0, 0);
1007	}
1008	BOOST_UBLAS_INLINE
1009	const_iterator1 end1 () const {
1010	return find1 (0, size1 (), 0);
1011	}
1012
1013	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1014	class const_iterator2:
1015	public container_const_reference<matrix_unary1>,
1016	public iterator_base_traits<typename E::const_iterator2::iterator_category>::template
1017	iterator_base<const_iterator2, value_type>::type {
1018	public:
1019	typedef typename E::const_iterator2::iterator_category iterator_category;
1020	typedef typename matrix_unary1::difference_type difference_type;
1021	typedef typename matrix_unary1::value_type value_type;
1022	typedef typename matrix_unary1::const_reference reference;
1023	typedef typename matrix_unary1::const_pointer pointer;
1024
1025	typedef const_iterator1 dual_iterator_type;
1026	typedef const_reverse_iterator1 dual_reverse_iterator_type;
1027
1028	// Construction and destruction
1029	BOOST_UBLAS_INLINE
1030	const_iterator2 ():
1031	container_const_reference<self_type> (), it_ () {}
1032	BOOST_UBLAS_INLINE
1033	const_iterator2 (const self_type &mu, const const_subiterator2_type &it):
1034	container_const_reference<self_type> (mu), it_ (it) {}
1035
1036	// Arithmetic
1037	BOOST_UBLAS_INLINE
1038	const_iterator2 &operator ++ () {
1039	++ it_;
1040	return *this;
1041	}
1042	BOOST_UBLAS_INLINE
1043	const_iterator2 &operator -- () {
1044	-- it_;
1045	return *this;
1046	}
1047	BOOST_UBLAS_INLINE
1048	const_iterator2 &operator += (difference_type n) {
1049	it_ += n;
1050	return *this;
1051	}
1052	BOOST_UBLAS_INLINE
1053	const_iterator2 &operator -= (difference_type n) {
1054	it_ -= n;
1055	return *this;
1056	}
1057	BOOST_UBLAS_INLINE
1058	difference_type operator - (const const_iterator2 &it) const {
1059	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1060	return it_ - it.it_;
1061	}
1062
1063	// Dereference
1064	BOOST_UBLAS_INLINE
1065	const_reference operator * () const {
1066	return functor_type::apply (*it_);
1067	}
1068
1069	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1070	BOOST_UBLAS_INLINE
1071	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1072	typename self_type::
1073	#endif
1074	const_iterator1 begin () const {
1075	return (*this) ().find1 (1, 0, index2 ());
1076	}
1077	BOOST_UBLAS_INLINE
1078	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1079	typename self_type::
1080	#endif
1081	const_iterator1 end () const {
1082	return (this) ().find1 (1, (this) ().size1 (), index2 ());
1083	}
1084	BOOST_UBLAS_INLINE
1085	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1086	typename self_type::
1087	#endif
1088	const_reverse_iterator1 rbegin () const {
1089	return const_reverse_iterator1 (end ());
1090	}
1091	BOOST_UBLAS_INLINE
1092	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1093	typename self_type::
1094	#endif
1095	const_reverse_iterator1 rend () const {
1096	return const_reverse_iterator1 (begin ());
1097	}
1098	#endif
1099
1100	// Indices
1101	BOOST_UBLAS_INLINE
1102	size_type index1 () const {
1103	return it_.index1 ();
1104	}
1105	BOOST_UBLAS_INLINE
1106	size_type index2 () const {
1107	return it_.index2 ();
1108	}
1109
1110	// Assignment
1111	BOOST_UBLAS_INLINE
1112	const_iterator2 &operator = (const const_iterator2 &it) {
1113	container_const_reference<self_type>::assign (&it ());
1114	it_ = it.it_;
1115	return *this;
1116	}
1117
1118	// Comparison
1119	BOOST_UBLAS_INLINE
1120	bool operator == (const const_iterator2 &it) const {
1121	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1122	return it_ == it.it_;
1123	}
1124	BOOST_UBLAS_INLINE
1125	bool operator < (const const_iterator2 &it) const {
1126	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1127	return it_ < it.it_;
1128	}
1129
1130	private:
1131	const_subiterator2_type it_;
1132	};
1133	#endif
1134
1135	BOOST_UBLAS_INLINE
1136	const_iterator2 begin2 () const {
1137	return find2 (0, 0, 0);
1138	}
1139	BOOST_UBLAS_INLINE
1140	const_iterator2 end2 () const {
1141	return find2 (0, 0, size2 ());
1142	}
1143
1144	// Reverse iterators
1145
1146	BOOST_UBLAS_INLINE
1147	const_reverse_iterator1 rbegin1 () const {
1148	return const_reverse_iterator1 (end1 ());
1149	}
1150	BOOST_UBLAS_INLINE
1151	const_reverse_iterator1 rend1 () const {
1152	return const_reverse_iterator1 (begin1 ());
1153	}
1154
1155	BOOST_UBLAS_INLINE
1156	const_reverse_iterator2 rbegin2 () const {
1157	return const_reverse_iterator2 (end2 ());
1158	}
1159	BOOST_UBLAS_INLINE
1160	const_reverse_iterator2 rend2 () const {
1161	return const_reverse_iterator2 (begin2 ());
1162	}
1163
1164	private:
1165	expression_closure_type e_;
1166	};
1167
1168	template<class E, class F>
1169	struct matrix_unary1_traits {
1170	typedef matrix_unary1<E, F> expression_type;
1171	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1172	typedef expression_type result_type;
1173	#else
1174	typedef typename E::matrix_temporary_type result_type;
1175	#endif
1176	};
1177
1178	// (- m) [i] [j] = - m [i] [j]
1179	template<class E>
1180	BOOST_UBLAS_INLINE
1181	typename matrix_unary1_traits<E, scalar_negate<typename E::value_type> >::result_type
1182	operator - (const matrix_expression<E> &e) {
1183	typedef typename matrix_unary1_traits<E, scalar_negate<typename E::value_type> >::expression_type expression_type;
1184	return expression_type (e ());
1185	}
1186
1187	// (conj m) [i] [j] = conj (m [i] [j])
1188	template<class E>
1189	BOOST_UBLAS_INLINE
1190	typename matrix_unary1_traits<E, scalar_conj<typename E::value_type> >::result_type
1191	conj (const matrix_expression<E> &e) {
1192	typedef typename matrix_unary1_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
1193	return expression_type (e ());
1194	}
1195
1196	// (real m) [i] [j] = real (m [i] [j])
1197	template<class E>
1198	BOOST_UBLAS_INLINE
1199	typename matrix_unary1_traits<E, scalar_real<typename E::value_type> >::result_type
1200	real (const matrix_expression<E> &e) {
1201	typedef typename matrix_unary1_traits<E, scalar_real<typename E::value_type> >::expression_type expression_type;
1202	return expression_type (e ());
1203	}
1204
1205	// (imag m) [i] [j] = imag (m [i] [j])
1206	template<class E>
1207	BOOST_UBLAS_INLINE
1208	typename matrix_unary1_traits<E, scalar_imag<typename E::value_type> >::result_type
1209	imag (const matrix_expression<E> &e) {
1210	typedef typename matrix_unary1_traits<E, scalar_imag<typename E::value_type> >::expression_type expression_type;
1211	return expression_type (e ());
1212	}
1213
1214	template<class E, class F>
1215	class matrix_unary2:
1216	public matrix_expression<matrix_unary2<E, F> > {
1217
1218	typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<typename E::value_type> >,
1219	E,
1220	const E>::type expression_type;
1221	typedef F functor_type;
1222	public:
1223	typedef typename boost::mpl::if_<boost::is_const<expression_type>,
1224	typename E::const_closure_type,
1225	typename E::closure_type>::type expression_closure_type;
1226	private:
1227	typedef matrix_unary2<E, F> self_type;
1228	public:
1229	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1230	using matrix_expression<self_type>::operator ();
1231	#endif
1232	typedef typename E::size_type size_type;
1233	typedef typename E::difference_type difference_type;
1234	typedef typename F::result_type value_type;
1235	typedef value_type const_reference;
1236	typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<value_type> >,
1237	typename E::reference,
1238	value_type>::type reference;
1239
1240	typedef const self_type const_closure_type;
1241	typedef self_type closure_type;
1242	// typedef typename E::orientation_category orientation_category;
1243	typedef typename boost::mpl::if_<boost::is_same<typename E::orientation_category,
1244	row_major_tag>,
1245	column_major_tag,
1246	typename boost::mpl::if_<boost::is_same<typename E::orientation_category,
1247	column_major_tag>,
1248	row_major_tag,
1249	typename E::orientation_category>::type>::type orientation_category;
1250	// typedef unknown_storage_tag storage_category;
1251	typedef typename E::storage_category storage_category;
1252
1253	// Construction and destruction
1254	BOOST_UBLAS_INLINE
1255	// ISSUE may be used as mutable expression.
1256	// matrix_unary2 (const expression_type &e):
1257	explicit matrix_unary2 (expression_type &e):
1258	e_ (e) {}
1259
1260	// Accessors
1261	BOOST_UBLAS_INLINE
1262	size_type size1 () const {
1263	return e_.size2 ();
1264	}
1265	BOOST_UBLAS_INLINE
1266	size_type size2 () const {
1267	return e_.size1 ();
1268	}
1269
1270	public:
1271	// Expression accessors
1272	BOOST_UBLAS_INLINE
1273	const expression_closure_type &expression () const {
1274	return e_;
1275	}
1276
1277	public:
1278	// Element access
1279	BOOST_UBLAS_INLINE
1280	const_reference operator () (size_type i, size_type j) const {
1281	return functor_type::apply (e_ (j, i));
1282	}
1283	BOOST_UBLAS_INLINE
1284	reference operator () (size_type i, size_type j) {
1285	BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
1286	return e_ (j, i);
1287	}
1288
1289	// Closure comparison
1290	BOOST_UBLAS_INLINE
1291	bool same_closure (const matrix_unary2 &mu2) const {
1292	return (*this).expression ().same_closure (mu2.expression ());
1293	}
1294
1295	// Iterator types
1296	private:
1297	typedef typename E::const_iterator1 const_subiterator2_type;
1298	typedef typename E::const_iterator2 const_subiterator1_type;
1299	typedef const value_type *const_pointer;
1300
1301	public:
1302	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1303	typedef indexed_const_iterator1<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator1;
1304	typedef const_iterator1 iterator1;
1305	typedef indexed_const_iterator2<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator2;
1306	typedef const_iterator2 iterator2;
1307	#else
1308	class const_iterator1;
1309	typedef const_iterator1 iterator1;
1310	class const_iterator2;
1311	typedef const_iterator2 iterator2;
1312	#endif
1313	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1314	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1315
1316	// Element lookup
1317	BOOST_UBLAS_INLINE
1318	const_iterator1 find1 (int rank, size_type i, size_type j) const {
1319	const_subiterator1_type it1 (e_.find2 (rank, j, i));
1320	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1321	return const_iterator1 (*this, it1.index2 (), it1.index1 ());
1322	#else
1323	return const_iterator1 (*this, it1);
1324	#endif
1325	}
1326	BOOST_UBLAS_INLINE
1327	const_iterator2 find2 (int rank, size_type i, size_type j) const {
1328	const_subiterator2_type it2 (e_.find1 (rank, j, i));
1329	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1330	return const_iterator2 (*this, it2.index2 (), it2.index1 ());
1331	#else
1332	return const_iterator2 (*this, it2);
1333	#endif
1334	}
1335
1336	// Iterators enhance the iterators of the referenced expression
1337	// with the unary functor.
1338
1339	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1340	class const_iterator1:
1341	public container_const_reference<matrix_unary2>,
1342	public iterator_base_traits<typename E::const_iterator2::iterator_category>::template
1343	iterator_base<const_iterator1, value_type>::type {
1344	public:
1345	typedef typename E::const_iterator2::iterator_category iterator_category;
1346	typedef typename matrix_unary2::difference_type difference_type;
1347	typedef typename matrix_unary2::value_type value_type;
1348	typedef typename matrix_unary2::const_reference reference;
1349	typedef typename matrix_unary2::const_pointer pointer;
1350
1351	typedef const_iterator2 dual_iterator_type;
1352	typedef const_reverse_iterator2 dual_reverse_iterator_type;
1353
1354	// Construction and destruction
1355	BOOST_UBLAS_INLINE
1356	const_iterator1 ():
1357	container_const_reference<self_type> (), it_ () {}
1358	BOOST_UBLAS_INLINE
1359	const_iterator1 (const self_type &mu, const const_subiterator1_type &it):
1360	container_const_reference<self_type> (mu), it_ (it) {}
1361
1362	// Arithmetic
1363	BOOST_UBLAS_INLINE
1364	const_iterator1 &operator ++ () {
1365	++ it_;
1366	return *this;
1367	}
1368	BOOST_UBLAS_INLINE
1369	const_iterator1 &operator -- () {
1370	-- it_;
1371	return *this;
1372	}
1373	BOOST_UBLAS_INLINE
1374	const_iterator1 &operator += (difference_type n) {
1375	it_ += n;
1376	return *this;
1377	}
1378	BOOST_UBLAS_INLINE
1379	const_iterator1 &operator -= (difference_type n) {
1380	it_ -= n;
1381	return *this;
1382	}
1383	BOOST_UBLAS_INLINE
1384	difference_type operator - (const const_iterator1 &it) const {
1385	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1386	return it_ - it.it_;
1387	}
1388
1389	// Dereference
1390	BOOST_UBLAS_INLINE
1391	const_reference operator * () const {
1392	return functor_type::apply (*it_);
1393	}
1394
1395	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1396	BOOST_UBLAS_INLINE
1397	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1398	typename self_type::
1399	#endif
1400	const_iterator2 begin () const {
1401	return (*this) ().find2 (1, index1 (), 0);
1402	}
1403	BOOST_UBLAS_INLINE
1404	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1405	typename self_type::
1406	#endif
1407	const_iterator2 end () const {
1408	return (this) ().find2 (1, index1 (), (this) ().size2 ());
1409	}
1410	BOOST_UBLAS_INLINE
1411	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1412	typename self_type::
1413	#endif
1414	const_reverse_iterator2 rbegin () const {
1415	return const_reverse_iterator2 (end ());
1416	}
1417	BOOST_UBLAS_INLINE
1418	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1419	typename self_type::
1420	#endif
1421	const_reverse_iterator2 rend () const {
1422	return const_reverse_iterator2 (begin ());
1423	}
1424	#endif
1425
1426	// Indices
1427	BOOST_UBLAS_INLINE
1428	size_type index1 () const {
1429	return it_.index2 ();
1430	}
1431	BOOST_UBLAS_INLINE
1432	size_type index2 () const {
1433	return it_.index1 ();
1434	}
1435
1436	// Assignment
1437	BOOST_UBLAS_INLINE
1438	const_iterator1 &operator = (const const_iterator1 &it) {
1439	container_const_reference<self_type>::assign (&it ());
1440	it_ = it.it_;
1441	return *this;
1442	}
1443
1444	// Comparison
1445	BOOST_UBLAS_INLINE
1446	bool operator == (const const_iterator1 &it) const {
1447	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1448	return it_ == it.it_;
1449	}
1450	BOOST_UBLAS_INLINE
1451	bool operator < (const const_iterator1 &it) const {
1452	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1453	return it_ < it.it_;
1454	}
1455
1456	private:
1457	const_subiterator1_type it_;
1458	};
1459	#endif
1460
1461	BOOST_UBLAS_INLINE
1462	const_iterator1 begin1 () const {
1463	return find1 (0, 0, 0);
1464	}
1465	BOOST_UBLAS_INLINE
1466	const_iterator1 end1 () const {
1467	return find1 (0, size1 (), 0);
1468	}
1469
1470	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1471	class const_iterator2:
1472	public container_const_reference<matrix_unary2>,
1473	public iterator_base_traits<typename E::const_iterator1::iterator_category>::template
1474	iterator_base<const_iterator2, value_type>::type {
1475	public:
1476	typedef typename E::const_iterator1::iterator_category iterator_category;
1477	typedef typename matrix_unary2::difference_type difference_type;
1478	typedef typename matrix_unary2::value_type value_type;
1479	typedef typename matrix_unary2::const_reference reference;
1480	typedef typename matrix_unary2::const_pointer pointer;
1481
1482	typedef const_iterator1 dual_iterator_type;
1483	typedef const_reverse_iterator1 dual_reverse_iterator_type;
1484
1485	// Construction and destruction
1486	BOOST_UBLAS_INLINE
1487	const_iterator2 ():
1488	container_const_reference<self_type> (), it_ () {}
1489	BOOST_UBLAS_INLINE
1490	const_iterator2 (const self_type &mu, const const_subiterator2_type &it):
1491	container_const_reference<self_type> (mu), it_ (it) {}
1492
1493	// Arithmetic
1494	BOOST_UBLAS_INLINE
1495	const_iterator2 &operator ++ () {
1496	++ it_;
1497	return *this;
1498	}
1499	BOOST_UBLAS_INLINE
1500	const_iterator2 &operator -- () {
1501	-- it_;
1502	return *this;
1503	}
1504	BOOST_UBLAS_INLINE
1505	const_iterator2 &operator += (difference_type n) {
1506	it_ += n;
1507	return *this;
1508	}
1509	BOOST_UBLAS_INLINE
1510	const_iterator2 &operator -= (difference_type n) {
1511	it_ -= n;
1512	return *this;
1513	}
1514	BOOST_UBLAS_INLINE
1515	difference_type operator - (const const_iterator2 &it) const {
1516	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1517	return it_ - it.it_;
1518	}
1519
1520	// Dereference
1521	BOOST_UBLAS_INLINE
1522	const_reference operator * () const {
1523	return functor_type::apply (*it_);
1524	}
1525
1526	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1527	BOOST_UBLAS_INLINE
1528	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1529	typename self_type::
1530	#endif
1531	const_iterator1 begin () const {
1532	return (*this) ().find1 (1, 0, index2 ());
1533	}
1534	BOOST_UBLAS_INLINE
1535	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1536	typename self_type::
1537	#endif
1538	const_iterator1 end () const {
1539	return (this) ().find1 (1, (this) ().size1 (), index2 ());
1540	}
1541	BOOST_UBLAS_INLINE
1542	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1543	typename self_type::
1544	#endif
1545	const_reverse_iterator1 rbegin () const {
1546	return const_reverse_iterator1 (end ());
1547	}
1548	BOOST_UBLAS_INLINE
1549	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1550	typename self_type::
1551	#endif
1552	const_reverse_iterator1 rend () const {
1553	return const_reverse_iterator1 (begin ());
1554	}
1555	#endif
1556
1557	// Indices
1558	BOOST_UBLAS_INLINE
1559	size_type index1 () const {
1560	return it_.index2 ();
1561	}
1562	BOOST_UBLAS_INLINE
1563	size_type index2 () const {
1564	return it_.index1 ();
1565	}
1566
1567	// Assignment
1568	BOOST_UBLAS_INLINE
1569	const_iterator2 &operator = (const const_iterator2 &it) {
1570	container_const_reference<self_type>::assign (&it ());
1571	it_ = it.it_;
1572	return *this;
1573	}
1574
1575	// Comparison
1576	BOOST_UBLAS_INLINE
1577	bool operator == (const const_iterator2 &it) const {
1578	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1579	return it_ == it.it_;
1580	}
1581	BOOST_UBLAS_INLINE
1582	bool operator < (const const_iterator2 &it) const {
1583	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1584	return it_ < it.it_;
1585	}
1586
1587	private:
1588	const_subiterator2_type it_;
1589	};
1590	#endif
1591
1592	BOOST_UBLAS_INLINE
1593	const_iterator2 begin2 () const {
1594	return find2 (0, 0, 0);
1595	}
1596	BOOST_UBLAS_INLINE
1597	const_iterator2 end2 () const {
1598	return find2 (0, 0, size2 ());
1599	}
1600
1601	// Reverse iterators
1602
1603	BOOST_UBLAS_INLINE
1604	const_reverse_iterator1 rbegin1 () const {
1605	return const_reverse_iterator1 (end1 ());
1606	}
1607	BOOST_UBLAS_INLINE
1608	const_reverse_iterator1 rend1 () const {
1609	return const_reverse_iterator1 (begin1 ());
1610	}
1611
1612	BOOST_UBLAS_INLINE
1613	const_reverse_iterator2 rbegin2 () const {
1614	return const_reverse_iterator2 (end2 ());
1615	}
1616	BOOST_UBLAS_INLINE
1617	const_reverse_iterator2 rend2 () const {
1618	return const_reverse_iterator2 (begin2 ());
1619	}
1620
1621	private:
1622	expression_closure_type e_;
1623	};
1624
1625	template<class E, class F>
1626	struct matrix_unary2_traits {
1627	typedef matrix_unary2<E, F> expression_type;
1628	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1629	typedef expression_type result_type;
1630	#else
1631	typedef typename E::matrix_temporary_type result_type;
1632	#endif
1633	};
1634
1635	// (trans m) [i] [j] = m [j] [i]
1636	template<class E>
1637	BOOST_UBLAS_INLINE
1638	typename matrix_unary2_traits<const E, scalar_identity<typename E::value_type> >::result_type
1639	trans (const matrix_expression<E> &e) {
1640	typedef typename matrix_unary2_traits<const E, scalar_identity<typename E::value_type> >::expression_type expression_type;
1641	return expression_type (e ());
1642	}
1643	template<class E>
1644	BOOST_UBLAS_INLINE
1645	typename matrix_unary2_traits<E, scalar_identity<typename E::value_type> >::result_type
1646	trans (matrix_expression<E> &e) {
1647	typedef typename matrix_unary2_traits<E, scalar_identity<typename E::value_type> >::expression_type expression_type;
1648	return expression_type (e ());
1649	}
1650
1651	// (herm m) [i] [j] = conj (m [j] [i])
1652	template<class E>
1653	BOOST_UBLAS_INLINE
1654	typename matrix_unary2_traits<E, scalar_conj<typename E::value_type> >::result_type
1655	herm (const matrix_expression<E> &e) {
1656	typedef typename matrix_unary2_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
1657	return expression_type (e ());
1658	}
1659
1660	template<class E1, class E2, class F>
1661	class matrix_binary:
1662	public matrix_expression<matrix_binary<E1, E2, F> > {
1663
1664	typedef E1 expression1_type;
1665	typedef E2 expression2_type;
1666	typedef F functor_type;
1667	public:
1668	typedef typename E1::const_closure_type expression1_closure_type;
1669	typedef typename E2::const_closure_type expression2_closure_type;
1670	private:
1671	typedef matrix_binary<E1, E2, F> self_type;
1672	public:
1673	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1674	using matrix_expression<self_type>::operator ();
1675	#endif
1676	typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
1677	typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
1678	typedef typename F::result_type value_type;
1679	typedef value_type const_reference;
1680	typedef const_reference reference;
1681	typedef const self_type const_closure_type;
1682	typedef const_closure_type closure_type;
1683	typedef unknown_orientation_tag orientation_category;
1684	typedef unknown_storage_tag storage_category;
1685
1686	// Construction and destruction
1687	BOOST_UBLAS_INLINE
1688	matrix_binary (const E1 &e1, const E2 &e2):
1689	e1_ (e1), e2_ (e2) {}
1690
1691	// Accessors
1692	BOOST_UBLAS_INLINE
1693	size_type size1 () const {
1694	return BOOST_UBLAS_SAME (e1_.size1 (), e2_.size1 ());
1695	}
1696	BOOST_UBLAS_INLINE
1697	size_type size2 () const {
1698	return BOOST_UBLAS_SAME (e1_.size2 (), e2_.size2 ());
1699	}
1700
1701	public:
1702	// Expression accessors
1703	BOOST_UBLAS_INLINE
1704	const expression1_closure_type &expression1 () const {
1705	return e1_;
1706	}
1707	BOOST_UBLAS_INLINE
1708	const expression2_closure_type &expression2 () const {
1709	return e2_;
1710	}
1711
1712	public:
1713	// Element access
1714	BOOST_UBLAS_INLINE
1715	const_reference operator () (size_type i, size_type j) const {
1716	return functor_type::apply (e1_ (i, j), e2_ (i, j));
1717	}
1718
1719	// Closure comparison
1720	BOOST_UBLAS_INLINE
1721	bool same_closure (const matrix_binary &mb) const {
1722	return (*this).expression1 ().same_closure (mb.expression1 ()) &&
1723	(*this).expression2 ().same_closure (mb.expression2 ());
1724	}
1725
1726	// Iterator types
1727	private:
1728	typedef typename E1::const_iterator1 const_iterator11_type;
1729	typedef typename E1::const_iterator2 const_iterator12_type;
1730	typedef typename E2::const_iterator1 const_iterator21_type;
1731	typedef typename E2::const_iterator2 const_iterator22_type;
1732	typedef const value_type *const_pointer;
1733
1734	public:
1735	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1736	typedef typename iterator_restrict_traits<typename const_iterator11_type::iterator_category,
1737	typename const_iterator21_type::iterator_category>::iterator_category iterator_category1;
1738	typedef indexed_const_iterator1<const_closure_type, iterator_category1> const_iterator1;
1739	typedef const_iterator1 iterator1;
1740	typedef typename iterator_restrict_traits<typename const_iterator12_type::iterator_category,
1741	typename const_iterator22_type::iterator_category>::iterator_category iterator_category2;
1742	typedef indexed_const_iterator2<const_closure_type, iterator_category2> const_iterator2;
1743	typedef const_iterator2 iterator2;
1744	#else
1745	class const_iterator1;
1746	typedef const_iterator1 iterator1;
1747	class const_iterator2;
1748	typedef const_iterator2 iterator2;
1749	#endif
1750	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1751	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1752
1753	// Element lookup
1754	BOOST_UBLAS_INLINE
1755	const_iterator1 find1 (int rank, size_type i, size_type j) const {
1756	const_iterator11_type it11 (e1_.find1 (rank, i, j));
1757	const_iterator11_type it11_end (e1_.find1 (rank, size1 (), j));
1758	const_iterator21_type it21 (e2_.find1 (rank, i, j));
1759	const_iterator21_type it21_end (e2_.find1 (rank, size1 (), j));
1760	BOOST_UBLAS_CHECK (rank == 0 \|\| it11 == it11_end \|\| it11.index2 () == j, internal_logic ())
1761	BOOST_UBLAS_CHECK (rank == 0 \|\| it21 == it21_end \|\| it21.index2 () == j, internal_logic ())
1762	i = (std::min) (it11 != it11_end ? it11.index1 () : size1 (),
1763	it21 != it21_end ? it21.index1 () : size1 ());
1764	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1765	return const_iterator1 (*this, i, j);
1766	#else
1767	return const_iterator1 (*this, i, j, it11, it11_end, it21, it21_end);
1768	#endif
1769	}
1770	BOOST_UBLAS_INLINE
1771	const_iterator2 find2 (int rank, size_type i, size_type j) const {
1772	const_iterator12_type it12 (e1_.find2 (rank, i, j));
1773	const_iterator12_type it12_end (e1_.find2 (rank, i, size2 ()));
1774	const_iterator22_type it22 (e2_.find2 (rank, i, j));
1775	const_iterator22_type it22_end (e2_.find2 (rank, i, size2 ()));
1776	BOOST_UBLAS_CHECK (rank == 0 \|\| it12 == it12_end \|\| it12.index1 () == i, internal_logic ())
1777	BOOST_UBLAS_CHECK (rank == 0 \|\| it22 == it22_end \|\| it22.index1 () == i, internal_logic ())
1778	j = (std::min) (it12 != it12_end ? it12.index2 () : size2 (),
1779	it22 != it22_end ? it22.index2 () : size2 ());
1780	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1781	return const_iterator2 (*this, i, j);
1782	#else
1783	return const_iterator2 (*this, i, j, it12, it12_end, it22, it22_end);
1784	#endif
1785	}
1786
1787	// Iterators enhance the iterators of the referenced expression
1788	// with the binary functor.
1789
1790	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1791	class const_iterator1:
1792	public container_const_reference<matrix_binary>,
1793	public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
1794	typename E2::const_iterator1::iterator_category>::iterator_category>::template
1795	iterator_base<const_iterator1, value_type>::type {
1796	public:
1797	typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
1798	typename E2::const_iterator1::iterator_category>::iterator_category iterator_category;
1799	typedef typename matrix_binary::difference_type difference_type;
1800	typedef typename matrix_binary::value_type value_type;
1801	typedef typename matrix_binary::const_reference reference;
1802	typedef typename matrix_binary::const_pointer pointer;
1803
1804	typedef const_iterator2 dual_iterator_type;
1805	typedef const_reverse_iterator2 dual_reverse_iterator_type;
1806
1807	// Construction and destruction
1808	BOOST_UBLAS_INLINE
1809	const_iterator1 ():
1810	container_const_reference<self_type> (), i_ (), j_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
1811	BOOST_UBLAS_INLINE
1812	const_iterator1 (const self_type &mb, size_type i, size_type j,
1813	const const_iterator11_type &it1, const const_iterator11_type &it1_end,
1814	const const_iterator21_type &it2, const const_iterator21_type &it2_end):
1815	container_const_reference<self_type> (mb), i_ (i), j_ (j), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
1816
1817	// Dense specializations
1818	BOOST_UBLAS_INLINE
1819	void increment (dense_random_access_iterator_tag) {
1820	++ i_, ++ it1_, ++ it2_;
1821	}
1822	BOOST_UBLAS_INLINE
1823	void decrement (dense_random_access_iterator_tag) {
1824	-- i_, -- it1_, -- it2_;
1825	}
1826	BOOST_UBLAS_INLINE
1827	value_type dereference (dense_random_access_iterator_tag) const {
1828	return functor_type::apply (it1_, it2_);
1829	}
1830
1831	// Packed specializations
1832	BOOST_UBLAS_INLINE
1833	void increment (packed_random_access_iterator_tag) {
1834	if (it1_ != it1_end_)
1835	if (it1_.index1 () <= i_)
1836	++ it1_;
1837	if (it2_ != it2_end_)
1838	if (it2_.index1 () <= i_)
1839	++ it2_;
1840	++ i_;
1841	}
1842	BOOST_UBLAS_INLINE
1843	void decrement (packed_random_access_iterator_tag) {
1844	if (it1_ != it1_end_)
1845	if (i_ <= it1_.index1 ())
1846	-- it1_;
1847	if (it2_ != it2_end_)
1848	if (i_ <= it2_.index1 ())
1849	-- it2_;
1850	-- i_;
1851	}
1852	BOOST_UBLAS_INLINE
1853	value_type dereference (packed_random_access_iterator_tag) const {
1854	value_type t1 = value_type/zero/();
1855	if (it1_ != it1_end_) {
1856	BOOST_UBLAS_CHECK (it1_.index2 () == j_, internal_logic ());
1857	if (it1_.index1 () == i_)
1858	t1 = *it1_;
1859	}
1860	value_type t2 = value_type/zero/();
1861	if (it2_ != it2_end_) {
1862	BOOST_UBLAS_CHECK (it2_.index2 () == j_, internal_logic ());
1863	if (it2_.index1 () == i_)
1864	t2 = *it2_;
1865	}
1866	return functor_type::apply (t1, t2);
1867	}
1868
1869	// Sparse specializations
1870	BOOST_UBLAS_INLINE
1871	void increment (sparse_bidirectional_iterator_tag) {
1872	size_type index1 = (*this) ().size1 ();
1873	if (it1_ != it1_end_) {
1874	if (it1_.index1 () <= i_)
1875	++ it1_;
1876	if (it1_ != it1_end_)
1877	index1 = it1_.index1 ();
1878	}
1879	size_type index2 = (*this) ().size1 ();
1880	if (it2_ != it2_end_)
1881	if (it2_.index1 () <= i_)
1882	++ it2_;
1883	if (it2_ != it2_end_) {
1884	index2 = it2_.index1 ();
1885	}
1886	i_ = (std::min) (index1, index2);
1887	}
1888	BOOST_UBLAS_INLINE
1889	void decrement (sparse_bidirectional_iterator_tag) {
1890	size_type index1 = (*this) ().size1 ();
1891	if (it1_ != it1_end_) {
1892	if (i_ <= it1_.index1 ())
1893	-- it1_;
1894	if (it1_ != it1_end_)
1895	index1 = it1_.index1 ();
1896	}
1897	size_type index2 = (*this) ().size1 ();
1898	if (it2_ != it2_end_) {
1899	if (i_ <= it2_.index1 ())
1900	-- it2_;
1901	if (it2_ != it2_end_)
1902	index2 = it2_.index1 ();
1903	}
1904	i_ = (std::max) (index1, index2);
1905	}
1906	BOOST_UBLAS_INLINE
1907	value_type dereference (sparse_bidirectional_iterator_tag) const {
1908	value_type t1 = value_type/zero/();
1909	if (it1_ != it1_end_) {
1910	BOOST_UBLAS_CHECK (it1_.index2 () == j_, internal_logic ());
1911	if (it1_.index1 () == i_)
1912	t1 = *it1_;
1913	}
1914	value_type t2 = value_type/zero/();
1915	if (it2_ != it2_end_) {
1916	BOOST_UBLAS_CHECK (it2_.index2 () == j_, internal_logic ());
1917	if (it2_.index1 () == i_)
1918	t2 = *it2_;
1919	}
1920	return functor_type::apply (t1, t2);
1921	}
1922
1923	// Arithmetic
1924	BOOST_UBLAS_INLINE
1925	const_iterator1 &operator ++ () {
1926	increment (iterator_category ());
1927	return *this;
1928	}
1929	BOOST_UBLAS_INLINE
1930	const_iterator1 &operator -- () {
1931	decrement (iterator_category ());
1932	return *this;
1933	}
1934	BOOST_UBLAS_INLINE
1935	const_iterator1 &operator += (difference_type n) {
1936	i_ += n, it1_ += n, it2_ += n;
1937	return *this;
1938	}
1939	BOOST_UBLAS_INLINE
1940	const_iterator1 &operator -= (difference_type n) {
1941	i_ -= n, it1_ -= n, it2_ -= n;
1942	return *this;
1943	}
1944	BOOST_UBLAS_INLINE
1945	difference_type operator - (const const_iterator1 &it) const {
1946	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1947	BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
1948	return index1 () - it.index1 ();
1949	}
1950
1951	// Dereference
1952	BOOST_UBLAS_INLINE
1953	const_reference operator * () const {
1954	return dereference (iterator_category ());
1955	}
1956
1957	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1958	BOOST_UBLAS_INLINE
1959	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1960	typename self_type::
1961	#endif
1962	const_iterator2 begin () const {
1963	return (*this) ().find2 (1, index1 (), 0);
1964	}
1965	BOOST_UBLAS_INLINE
1966	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1967	typename self_type::
1968	#endif
1969	const_iterator2 end () const {
1970	return (this) ().find2 (1, index1 (), (this) ().size2 ());
1971	}
1972	BOOST_UBLAS_INLINE
1973	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1974	typename self_type::
1975	#endif
1976	const_reverse_iterator2 rbegin () const {
1977	return const_reverse_iterator2 (end ());
1978	}
1979	BOOST_UBLAS_INLINE
1980	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1981	typename self_type::
1982	#endif
1983	const_reverse_iterator2 rend () const {
1984	return const_reverse_iterator2 (begin ());
1985	}
1986	#endif
1987
1988	// Indices
1989	BOOST_UBLAS_INLINE
1990	size_type index1 () const {
1991	return i_;
1992	}
1993	BOOST_UBLAS_INLINE
1994	size_type index2 () const {
1995	// if (it1_ != it1_end_ && it2_ != it2_end_)
1996	// return BOOST_UBLAS_SAME (it1_.index2 (), it2_.index2 ());
1997	// else
1998	return j_;
1999	}
2000
2001	// Assignment
2002	BOOST_UBLAS_INLINE
2003	const_iterator1 &operator = (const const_iterator1 &it) {
2004	container_const_reference<self_type>::assign (&it ());
2005	i_ = it.i_;
2006	j_ = it.j_;
2007	it1_ = it.it1_;
2008	it1_end_ = it.it1_end_;
2009	it2_ = it.it2_;
2010	it2_end_ = it.it2_end_;
2011	return *this;
2012	}
2013
2014	// Comparison
2015	BOOST_UBLAS_INLINE
2016	bool operator == (const const_iterator1 &it) const {
2017	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2018	BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
2019	return index1 () == it.index1 ();
2020	}
2021	BOOST_UBLAS_INLINE
2022	bool operator < (const const_iterator1 &it) const {
2023	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2024	BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
2025	return index1 () < it.index1 ();
2026	}
2027
2028	private:
2029	size_type i_;
2030	size_type j_;
2031	const_iterator11_type it1_;
2032	const_iterator11_type it1_end_;
2033	const_iterator21_type it2_;
2034	const_iterator21_type it2_end_;
2035	};
2036	#endif
2037
2038	BOOST_UBLAS_INLINE
2039	const_iterator1 begin1 () const {
2040	return find1 (0, 0, 0);
2041	}
2042	BOOST_UBLAS_INLINE
2043	const_iterator1 end1 () const {
2044	return find1 (0, size1 (), 0);
2045	}
2046
2047	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2048	class const_iterator2:
2049	public container_const_reference<matrix_binary>,
2050	public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator2::iterator_category,
2051	typename E2::const_iterator2::iterator_category>::iterator_category>::template
2052	iterator_base<const_iterator2, value_type>::type {
2053	public:
2054	typedef typename iterator_restrict_traits<typename E1::const_iterator2::iterator_category,
2055	typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
2056	typedef typename matrix_binary::difference_type difference_type;
2057	typedef typename matrix_binary::value_type value_type;
2058	typedef typename matrix_binary::const_reference reference;
2059	typedef typename matrix_binary::const_pointer pointer;
2060
2061	typedef const_iterator1 dual_iterator_type;
2062	typedef const_reverse_iterator1 dual_reverse_iterator_type;
2063
2064	// Construction and destruction
2065	BOOST_UBLAS_INLINE
2066	const_iterator2 ():
2067	container_const_reference<self_type> (), i_ (), j_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
2068	BOOST_UBLAS_INLINE
2069	const_iterator2 (const self_type &mb, size_type i, size_type j,
2070	const const_iterator12_type &it1, const const_iterator12_type &it1_end,
2071	const const_iterator22_type &it2, const const_iterator22_type &it2_end):
2072	container_const_reference<self_type> (mb), i_ (i), j_ (j), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
2073
2074	// Dense access specializations
2075	BOOST_UBLAS_INLINE
2076	void increment (dense_random_access_iterator_tag) {
2077	++ j_, ++ it1_, ++ it2_;
2078	}
2079	BOOST_UBLAS_INLINE
2080	void decrement (dense_random_access_iterator_tag) {
2081	-- j_, -- it1_, -- it2_;
2082	}
2083	BOOST_UBLAS_INLINE
2084	value_type dereference (dense_random_access_iterator_tag) const {
2085	return functor_type::apply (it1_, it2_);
2086	}
2087
2088	// Packed specializations
2089	BOOST_UBLAS_INLINE
2090	void increment (packed_random_access_iterator_tag) {
2091	if (it1_ != it1_end_)
2092	if (it1_.index2 () <= j_)
2093	++ it1_;
2094	if (it2_ != it2_end_)
2095	if (it2_.index2 () <= j_)
2096	++ it2_;
2097	++ j_;
2098	}
2099	BOOST_UBLAS_INLINE
2100	void decrement (packed_random_access_iterator_tag) {
2101	if (it1_ != it1_end_)
2102	if (j_ <= it1_.index2 ())
2103	-- it1_;
2104	if (it2_ != it2_end_)
2105	if (j_ <= it2_.index2 ())
2106	-- it2_;
2107	-- j_;
2108	}
2109	BOOST_UBLAS_INLINE
2110	value_type dereference (packed_random_access_iterator_tag) const {
2111	value_type t1 = value_type/zero/();
2112	if (it1_ != it1_end_) {
2113	BOOST_UBLAS_CHECK (it1_.index1 () == i_, internal_logic ());
2114	if (it1_.index2 () == j_)
2115	t1 = *it1_;
2116	}
2117	value_type t2 = value_type/zero/();
2118	if (it2_ != it2_end_) {
2119	BOOST_UBLAS_CHECK (it2_.index1 () == i_, internal_logic ());
2120	if (it2_.index2 () == j_)
2121	t2 = *it2_;
2122	}
2123	return functor_type::apply (t1, t2);
2124	}
2125
2126	// Sparse specializations
2127	BOOST_UBLAS_INLINE
2128	void increment (sparse_bidirectional_iterator_tag) {
2129	size_type index1 = (*this) ().size2 ();
2130	if (it1_ != it1_end_) {
2131	if (it1_.index2 () <= j_)
2132	++ it1_;
2133	if (it1_ != it1_end_)
2134	index1 = it1_.index2 ();
2135	}
2136	size_type index2 = (*this) ().size2 ();
2137	if (it2_ != it2_end_) {
2138	if (it2_.index2 () <= j_)
2139	++ it2_;
2140	if (it2_ != it2_end_)
2141	index2 = it2_.index2 ();
2142	}
2143	j_ = (std::min) (index1, index2);
2144	}
2145	BOOST_UBLAS_INLINE
2146	void decrement (sparse_bidirectional_iterator_tag) {
2147	size_type index1 = (*this) ().size2 ();
2148	if (it1_ != it1_end_) {
2149	if (j_ <= it1_.index2 ())
2150	-- it1_;
2151	if (it1_ != it1_end_)
2152	index1 = it1_.index2 ();
2153	}
2154	size_type index2 = (*this) ().size2 ();
2155	if (it2_ != it2_end_) {
2156	if (j_ <= it2_.index2 ())
2157	-- it2_;
2158	if (it2_ != it2_end_)
2159	index2 = it2_.index2 ();
2160	}
2161	j_ = (std::max) (index1, index2);
2162	}
2163	BOOST_UBLAS_INLINE
2164	value_type dereference (sparse_bidirectional_iterator_tag) const {
2165	value_type t1 = value_type/zero/();
2166	if (it1_ != it1_end_) {
2167	BOOST_UBLAS_CHECK (it1_.index1 () == i_, internal_logic ());
2168	if (it1_.index2 () == j_)
2169	t1 = *it1_;
2170	}
2171	value_type t2 = value_type/zero/();
2172	if (it2_ != it2_end_) {
2173	BOOST_UBLAS_CHECK (it2_.index1 () == i_, internal_logic ());
2174	if (it2_.index2 () == j_)
2175	t2 = *it2_;
2176	}
2177	return functor_type::apply (t1, t2);
2178	}
2179
2180	// Arithmetic
2181	BOOST_UBLAS_INLINE
2182	const_iterator2 &operator ++ () {
2183	increment (iterator_category ());
2184	return *this;
2185	}
2186	BOOST_UBLAS_INLINE
2187	const_iterator2 &operator -- () {
2188	decrement (iterator_category ());
2189	return *this;
2190	}
2191	BOOST_UBLAS_INLINE
2192	const_iterator2 &operator += (difference_type n) {
2193	j_ += n, it1_ += n, it2_ += n;
2194	return *this;
2195	}
2196	BOOST_UBLAS_INLINE
2197	const_iterator2 &operator -= (difference_type n) {
2198	j_ -= n, it1_ -= n, it2_ -= n;
2199	return *this;
2200	}
2201	BOOST_UBLAS_INLINE
2202	difference_type operator - (const const_iterator2 &it) const {
2203	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2204	BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2205	return index2 () - it.index2 ();
2206	}
2207
2208	// Dereference
2209	BOOST_UBLAS_INLINE
2210	const_reference operator * () const {
2211	return dereference (iterator_category ());
2212	}
2213
2214	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2215	BOOST_UBLAS_INLINE
2216	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2217	typename self_type::
2218	#endif
2219	const_iterator1 begin () const {
2220	return (*this) ().find1 (1, 0, index2 ());
2221	}
2222	BOOST_UBLAS_INLINE
2223	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2224	typename self_type::
2225	#endif
2226	const_iterator1 end () const {
2227	return (this) ().find1 (1, (this) ().size1 (), index2 ());
2228	}
2229	BOOST_UBLAS_INLINE
2230	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2231	typename self_type::
2232	#endif
2233	const_reverse_iterator1 rbegin () const {
2234	return const_reverse_iterator1 (end ());
2235	}
2236	BOOST_UBLAS_INLINE
2237	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2238	typename self_type::
2239	#endif
2240	const_reverse_iterator1 rend () const {
2241	return const_reverse_iterator1 (begin ());
2242	}
2243	#endif
2244
2245	// Indices
2246	BOOST_UBLAS_INLINE
2247	size_type index1 () const {
2248	// if (it1_ != it1_end_ && it2_ != it2_end_)
2249	// return BOOST_UBLAS_SAME (it1_.index1 (), it2_.index1 ());
2250	// else
2251	return i_;
2252	}
2253	BOOST_UBLAS_INLINE
2254	size_type index2 () const {
2255	return j_;
2256	}
2257
2258	// Assignment
2259	BOOST_UBLAS_INLINE
2260	const_iterator2 &operator = (const const_iterator2 &it) {
2261	container_const_reference<self_type>::assign (&it ());
2262	i_ = it.i_;
2263	j_ = it.j_;
2264	it1_ = it.it1_;
2265	it1_end_ = it.it1_end_;
2266	it2_ = it.it2_;
2267	it2_end_ = it.it2_end_;
2268	return *this;
2269	}
2270
2271	// Comparison
2272	BOOST_UBLAS_INLINE
2273	bool operator == (const const_iterator2 &it) const {
2274	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2275	BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2276	return index2 () == it.index2 ();
2277	}
2278	BOOST_UBLAS_INLINE
2279	bool operator < (const const_iterator2 &it) const {
2280	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2281	BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2282	return index2 () < it.index2 ();
2283	}
2284
2285	private:
2286	size_type i_;
2287	size_type j_;
2288	const_iterator12_type it1_;
2289	const_iterator12_type it1_end_;
2290	const_iterator22_type it2_;
2291	const_iterator22_type it2_end_;
2292	};
2293	#endif
2294
2295	BOOST_UBLAS_INLINE
2296	const_iterator2 begin2 () const {
2297	return find2 (0, 0, 0);
2298	}
2299	BOOST_UBLAS_INLINE
2300	const_iterator2 end2 () const {
2301	return find2 (0, 0, size2 ());
2302	}
2303
2304	// Reverse iterators
2305
2306	BOOST_UBLAS_INLINE
2307	const_reverse_iterator1 rbegin1 () const {
2308	return const_reverse_iterator1 (end1 ());
2309	}
2310	BOOST_UBLAS_INLINE
2311	const_reverse_iterator1 rend1 () const {
2312	return const_reverse_iterator1 (begin1 ());
2313	}
2314
2315	BOOST_UBLAS_INLINE
2316	const_reverse_iterator2 rbegin2 () const {
2317	return const_reverse_iterator2 (end2 ());
2318	}
2319	BOOST_UBLAS_INLINE
2320	const_reverse_iterator2 rend2 () const {
2321	return const_reverse_iterator2 (begin2 ());
2322	}
2323
2324	private:
2325	expression1_closure_type e1_;
2326	expression2_closure_type e2_;
2327	};
2328
2329	template<class E1, class E2, class F>
2330	struct matrix_binary_traits {
2331	typedef matrix_binary<E1, E2, F> expression_type;
2332	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
2333	typedef expression_type result_type;
2334	#else
2335	typedef typename E1::matrix_temporary_type result_type;
2336	#endif
2337	};
2338
2339	// (m1 + m2) [i] [j] = m1 [i] [j] + m2 [i] [j]
2340	template<class E1, class E2>
2341	BOOST_UBLAS_INLINE
2342	typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
2343	typename E2::value_type> >::result_type
2344	operator + (const matrix_expression<E1> &e1,
2345	const matrix_expression<E2> &e2) {
2346	typedef typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
2347	typename E2::value_type> >::expression_type expression_type;
2348	return expression_type (e1 (), e2 ());
2349	}
2350
2351	// (m1 - m2) [i] [j] = m1 [i] [j] - m2 [i] [j]
2352	template<class E1, class E2>
2353	BOOST_UBLAS_INLINE
2354	typename matrix_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
2355	typename E2::value_type> >::result_type
2356	operator - (const matrix_expression<E1> &e1,
2357	const matrix_expression<E2> &e2) {
2358	typedef typename matrix_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
2359	typename E2::value_type> >::expression_type expression_type;
2360	return expression_type (e1 (), e2 ());
2361	}
2362
2363	// (m1 * m2) [i] [j] = m1 [i] [j] * m2 [i] [j]
2364	template<class E1, class E2>
2365	BOOST_UBLAS_INLINE
2366	typename matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
2367	typename E2::value_type> >::result_type
2368	element_prod (const matrix_expression<E1> &e1,
2369	const matrix_expression<E2> &e2) {
2370	typedef typename matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
2371	typename E2::value_type> >::expression_type expression_type;
2372	return expression_type (e1 (), e2 ());
2373	}
2374
2375	// (m1 / m2) [i] [j] = m1 [i] [j] / m2 [i] [j]
2376	template<class E1, class E2>
2377	BOOST_UBLAS_INLINE
2378	typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
2379	typename E2::value_type> >::result_type
2380	element_div (const matrix_expression<E1> &e1,
2381	const matrix_expression<E2> &e2) {
2382	typedef typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
2383	typename E2::value_type> >::expression_type expression_type;
2384	return expression_type (e1 (), e2 ());
2385	}
2386
2387	template<class E1, class E2, class F>
2388	class matrix_binary_scalar1:
2389	public matrix_expression<matrix_binary_scalar1<E1, E2, F> > {
2390
2391	typedef E1 expression1_type;
2392	typedef E2 expression2_type;
2393	typedef F functor_type;
2394	typedef const E1& expression1_closure_type;
2395	typedef typename E2::const_closure_type expression2_closure_type;
2396	typedef matrix_binary_scalar1<E1, E2, F> self_type;
2397	public:
2398	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
2399	using matrix_expression<self_type>::operator ();
2400	#endif
2401	typedef typename E2::size_type size_type;
2402	typedef typename E2::difference_type difference_type;
2403	typedef typename F::result_type value_type;
2404	typedef value_type const_reference;
2405	typedef const_reference reference;
2406	typedef const self_type const_closure_type;
2407	typedef const_closure_type closure_type;
2408	typedef typename E2::orientation_category orientation_category;
2409	typedef unknown_storage_tag storage_category;
2410
2411	// Construction and destruction
2412	BOOST_UBLAS_INLINE
2413	matrix_binary_scalar1 (const expression1_type &e1, const expression2_type &e2):
2414	e1_ (e1), e2_ (e2) {}
2415
2416	// Accessors
2417	BOOST_UBLAS_INLINE
2418	size_type size1 () const {
2419	return e2_.size1 ();
2420	}
2421	BOOST_UBLAS_INLINE
2422	size_type size2 () const {
2423	return e2_.size2 ();
2424	}
2425
2426	public:
2427	// Element access
2428	BOOST_UBLAS_INLINE
2429	const_reference operator () (size_type i, size_type j) const {
2430	return functor_type::apply (expression1_type (e1_), e2_ (i, j));
2431	}
2432
2433	// Closure comparison
2434	BOOST_UBLAS_INLINE
2435	bool same_closure (const matrix_binary_scalar1 &mbs1) const {
2436	return &e1_ == &(mbs1.e1_) &&
2437	(*this).e2_.same_closure (mbs1.e2_);
2438	}
2439
2440	// Iterator types
2441	private:
2442	typedef expression1_type const_subiterator1_type;
2443	typedef typename E2::const_iterator1 const_iterator21_type;
2444	typedef typename E2::const_iterator2 const_iterator22_type;
2445	typedef const value_type *const_pointer;
2446
2447	public:
2448	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2449	typedef indexed_const_iterator1<const_closure_type, typename const_iterator21_type::iterator_category> const_iterator1;
2450	typedef const_iterator1 iterator1;
2451	typedef indexed_const_iterator2<const_closure_type, typename const_iterator22_type::iterator_category> const_iterator2;
2452	typedef const_iterator2 iterator2;
2453	#else
2454	class const_iterator1;
2455	typedef const_iterator1 iterator1;
2456	class const_iterator2;
2457	typedef const_iterator2 iterator2;
2458	#endif
2459	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
2460	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
2461
2462	// Element lookup
2463	BOOST_UBLAS_INLINE
2464	const_iterator1 find1 (int rank, size_type i, size_type j) const {
2465	const_iterator21_type it21 (e2_.find1 (rank, i, j));
2466	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2467	return const_iterator1 (*this, it21.index1 (), it21.index2 ());
2468	#else
2469	return const_iterator1 (*this, const_subiterator1_type (e1_), it21);
2470	#endif
2471	}
2472	BOOST_UBLAS_INLINE
2473	const_iterator2 find2 (int rank, size_type i, size_type j) const {
2474	const_iterator22_type it22 (e2_.find2 (rank, i, j));
2475	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2476	return const_iterator2 (*this, it22.index1 (), it22.index2 ());
2477	#else
2478	return const_iterator2 (*this, const_subiterator1_type (e1_), it22);
2479	#endif
2480	}
2481
2482	// Iterators enhance the iterators of the referenced expression
2483	// with the binary functor.
2484
2485	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2486	class const_iterator1:
2487	public container_const_reference<matrix_binary_scalar1>,
2488	public iterator_base_traits<typename E2::const_iterator1::iterator_category>::template
2489	iterator_base<const_iterator1, value_type>::type {
2490	public:
2491	typedef typename E2::const_iterator1::iterator_category iterator_category;
2492	typedef typename matrix_binary_scalar1::difference_type difference_type;
2493	typedef typename matrix_binary_scalar1::value_type value_type;
2494	typedef typename matrix_binary_scalar1::const_reference reference;
2495	typedef typename matrix_binary_scalar1::const_pointer pointer;
2496
2497	typedef const_iterator2 dual_iterator_type;
2498	typedef const_reverse_iterator2 dual_reverse_iterator_type;
2499
2500	// Construction and destruction
2501	BOOST_UBLAS_INLINE
2502	const_iterator1 ():
2503	container_const_reference<self_type> (), it1_ (), it2_ () {}
2504	BOOST_UBLAS_INLINE
2505	const_iterator1 (const self_type &mbs, const const_subiterator1_type &it1, const const_iterator21_type &it2):
2506	container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
2507
2508	// Arithmetic
2509	BOOST_UBLAS_INLINE
2510	const_iterator1 &operator ++ () {
2511	++ it2_;
2512	return *this;
2513	}
2514	BOOST_UBLAS_INLINE
2515	const_iterator1 &operator -- () {
2516	-- it2_ ;
2517	return *this;
2518	}
2519	BOOST_UBLAS_INLINE
2520	const_iterator1 &operator += (difference_type n) {
2521	it2_ += n;
2522	return *this;
2523	}
2524	BOOST_UBLAS_INLINE
2525	const_iterator1 &operator -= (difference_type n) {
2526	it2_ -= n;
2527	return *this;
2528	}
2529	BOOST_UBLAS_INLINE
2530	difference_type operator - (const const_iterator1 &it) const {
2531	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2532	// FIXME we shouldn't compare floats
2533	// BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
2534	return it2_ - it.it2_;
2535	}
2536
2537	// Dereference
2538	BOOST_UBLAS_INLINE
2539	const_reference operator * () const {
2540	return functor_type::apply (it1_, *it2_);
2541	}
2542
2543	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2544	BOOST_UBLAS_INLINE
2545	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2546	typename self_type::
2547	#endif
2548	const_iterator2 begin () const {
2549	return (*this) ().find2 (1, index1 (), 0);
2550	}
2551	BOOST_UBLAS_INLINE
2552	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2553	typename self_type::
2554	#endif
2555	const_iterator2 end () const {
2556	return (this) ().find2 (1, index1 (), (this) ().size2 ());
2557	}
2558	BOOST_UBLAS_INLINE
2559	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2560	typename self_type::
2561	#endif
2562	const_reverse_iterator2 rbegin () const {
2563	return const_reverse_iterator2 (end ());
2564	}
2565	BOOST_UBLAS_INLINE
2566	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2567	typename self_type::
2568	#endif
2569	const_reverse_iterator2 rend () const {
2570	return const_reverse_iterator2 (begin ());
2571	}
2572	#endif
2573
2574	// Indices
2575	BOOST_UBLAS_INLINE
2576	size_type index1 () const {
2577	return it2_.index1 ();
2578	}
2579	BOOST_UBLAS_INLINE
2580	size_type index2 () const {
2581	return it2_.index2 ();
2582	}
2583
2584	// Assignment
2585	BOOST_UBLAS_INLINE
2586	const_iterator1 &operator = (const const_iterator1 &it) {
2587	container_const_reference<self_type>::assign (&it ());
2588	it1_ = it.it1_;
2589	it2_ = it.it2_;
2590	return *this;
2591	}
2592
2593	// Comparison
2594	BOOST_UBLAS_INLINE
2595	bool operator == (const const_iterator1 &it) const {
2596	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2597	// FIXME we shouldn't compare floats
2598	// BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
2599	return it2_ == it.it2_;
2600	}
2601	BOOST_UBLAS_INLINE
2602	bool operator < (const const_iterator1 &it) const {
2603	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2604	// FIXME we shouldn't compare floats
2605	// BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
2606	return it2_ < it.it2_;
2607	}
2608
2609	private:
2610	const_subiterator1_type it1_;
2611	const_iterator21_type it2_;
2612	};
2613	#endif
2614
2615	BOOST_UBLAS_INLINE
2616	const_iterator1 begin1 () const {
2617	return find1 (0, 0, 0);
2618	}
2619	BOOST_UBLAS_INLINE
2620	const_iterator1 end1 () const {
2621	return find1 (0, size1 (), 0);
2622	}
2623
2624	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2625	class const_iterator2:
2626	public container_const_reference<matrix_binary_scalar1>,
2627	public iterator_base_traits<typename E2::const_iterator2::iterator_category>::template
2628	iterator_base<const_iterator2, value_type>::type {
2629	public:
2630	typedef typename E2::const_iterator2::iterator_category iterator_category;
2631	typedef typename matrix_binary_scalar1::difference_type difference_type;
2632	typedef typename matrix_binary_scalar1::value_type value_type;
2633	typedef typename matrix_binary_scalar1::const_reference reference;
2634	typedef typename matrix_binary_scalar1::const_pointer pointer;
2635
2636	typedef const_iterator1 dual_iterator_type;
2637	typedef const_reverse_iterator1 dual_reverse_iterator_type;
2638
2639	// Construction and destruction
2640	BOOST_UBLAS_INLINE
2641	const_iterator2 ():
2642	container_const_reference<self_type> (), it1_ (), it2_ () {}
2643	BOOST_UBLAS_INLINE
2644	const_iterator2 (const self_type &mbs, const const_subiterator1_type &it1, const const_iterator22_type &it2):
2645	container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
2646
2647	// Arithmetic
2648	BOOST_UBLAS_INLINE
2649	const_iterator2 &operator ++ () {
2650	++ it2_;
2651	return *this;
2652	}
2653	BOOST_UBLAS_INLINE
2654	const_iterator2 &operator -- () {
2655	-- it2_;
2656	return *this;
2657	}
2658	BOOST_UBLAS_INLINE
2659	const_iterator2 &operator += (difference_type n) {
2660	it2_ += n;
2661	return *this;
2662	}
2663	BOOST_UBLAS_INLINE
2664	const_iterator2 &operator -= (difference_type n) {
2665	it2_ -= n;
2666	return *this;
2667	}
2668	BOOST_UBLAS_INLINE
2669	difference_type operator - (const const_iterator2 &it) const {
2670	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2671	// FIXME we shouldn't compare floats
2672	// BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
2673	return it2_ - it.it2_;
2674	}
2675
2676	// Dereference
2677	BOOST_UBLAS_INLINE
2678	const_reference operator * () const {
2679	return functor_type::apply (it1_, *it2_);
2680	}
2681
2682	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2683	BOOST_UBLAS_INLINE
2684	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2685	typename self_type::
2686	#endif
2687	const_iterator1 begin () const {
2688	return (*this) ().find1 (1, 0, index2 ());
2689	}
2690	BOOST_UBLAS_INLINE
2691	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2692	typename self_type::
2693	#endif
2694	const_iterator1 end () const {
2695	return (this) ().find1 (1, (this) ().size1 (), index2 ());
2696	}
2697	BOOST_UBLAS_INLINE
2698	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2699	typename self_type::
2700	#endif
2701	const_reverse_iterator1 rbegin () const {
2702	return const_reverse_iterator1 (end ());
2703	}
2704	BOOST_UBLAS_INLINE
2705	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2706	typename self_type::
2707	#endif
2708	const_reverse_iterator1 rend () const {
2709	return const_reverse_iterator1 (begin ());
2710	}
2711	#endif
2712
2713	// Indices
2714	BOOST_UBLAS_INLINE
2715	size_type index1 () const {
2716	return it2_.index1 ();
2717	}
2718	BOOST_UBLAS_INLINE
2719	size_type index2 () const {
2720	return it2_.index2 ();
2721	}
2722
2723	// Assignment
2724	BOOST_UBLAS_INLINE
2725	const_iterator2 &operator = (const const_iterator2 &it) {
2726	container_const_reference<self_type>::assign (&it ());
2727	it1_ = it.it1_;
2728	it2_ = it.it2_;
2729	return *this;
2730	}
2731
2732	// Comparison
2733	BOOST_UBLAS_INLINE
2734	bool operator == (const const_iterator2 &it) const {
2735	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2736	// FIXME we shouldn't compare floats
2737	// BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
2738	return it2_ == it.it2_;
2739	}
2740	BOOST_UBLAS_INLINE
2741	bool operator < (const const_iterator2 &it) const {
2742	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2743	// FIXME we shouldn't compare floats
2744	// BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
2745	return it2_ < it.it2_;
2746	}
2747
2748	private:
2749	const_subiterator1_type it1_;
2750	const_iterator22_type it2_;
2751	};
2752	#endif
2753
2754	BOOST_UBLAS_INLINE
2755	const_iterator2 begin2 () const {
2756	return find2 (0, 0, 0);
2757	}
2758	BOOST_UBLAS_INLINE
2759	const_iterator2 end2 () const {
2760	return find2 (0, 0, size2 ());
2761	}
2762
2763	// Reverse iterators
2764
2765	BOOST_UBLAS_INLINE
2766	const_reverse_iterator1 rbegin1 () const {
2767	return const_reverse_iterator1 (end1 ());
2768	}
2769	BOOST_UBLAS_INLINE
2770	const_reverse_iterator1 rend1 () const {
2771	return const_reverse_iterator1 (begin1 ());
2772	}
2773
2774	BOOST_UBLAS_INLINE
2775	const_reverse_iterator2 rbegin2 () const {
2776	return const_reverse_iterator2 (end2 ());
2777	}
2778	BOOST_UBLAS_INLINE
2779	const_reverse_iterator2 rend2 () const {
2780	return const_reverse_iterator2 (begin2 ());
2781	}
2782
2783	private:
2784	expression1_closure_type e1_;
2785	expression2_closure_type e2_;
2786	};
2787
2788	template<class E1, class E2, class F>
2789	struct matrix_binary_scalar1_traits {
2790	typedef matrix_binary_scalar1<E1, E2, F> expression_type; // allow E1 to be builtin type
2791	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
2792	typedef expression_type result_type;
2793	#else
2794	typedef typename E2::matrix_temporary_type result_type;
2795	#endif
2796	};
2797
2798	// (t * m) [i] [j] = t * m [i] [j]
2799	template<class T1, class E2>
2800	BOOST_UBLAS_INLINE
2801	typename matrix_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
2802	operator * (const T1 &e1,
2803	const matrix_expression<E2> &e2) {
2804	typedef typename matrix_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::expression_type expression_type;
2805	return expression_type (e1, e2 ());
2806	}
2807
2808
2809	template<class E1, class E2, class F>
2810	class matrix_binary_scalar2:
2811	public matrix_expression<matrix_binary_scalar2<E1, E2, F> > {
2812
2813	typedef E1 expression1_type;
2814	typedef E2 expression2_type;
2815	typedef F functor_type;
2816	public:
2817	typedef typename E1::const_closure_type expression1_closure_type;
2818	typedef const E2& expression2_closure_type;
2819	private:
2820	typedef matrix_binary_scalar2<E1, E2, F> self_type;
2821	public:
2822	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
2823	using matrix_expression<self_type>::operator ();
2824	#endif
2825	typedef typename E1::size_type size_type;
2826	typedef typename E1::difference_type difference_type;
2827	typedef typename F::result_type value_type;
2828	typedef value_type const_reference;
2829	typedef const_reference reference;
2830
2831	typedef const self_type const_closure_type;
2832	typedef const_closure_type closure_type;
2833	typedef typename E1::orientation_category orientation_category;
2834	typedef unknown_storage_tag storage_category;
2835
2836	// Construction and destruction
2837	BOOST_UBLAS_INLINE
2838	matrix_binary_scalar2 (const expression1_type &e1, const expression2_type &e2):
2839	e1_ (e1), e2_ (e2) {}
2840
2841	// Accessors
2842	BOOST_UBLAS_INLINE
2843	size_type size1 () const {
2844	return e1_.size1 ();
2845	}
2846	BOOST_UBLAS_INLINE
2847	size_type size2 () const {
2848	return e1_.size2 ();
2849	}
2850
2851	public:
2852	// Element access
2853	BOOST_UBLAS_INLINE
2854	const_reference operator () (size_type i, size_type j) const {
2855	return functor_type::apply (e1_ (i, j), expression2_type (e2_));
2856	}
2857
2858	// Closure comparison
2859	BOOST_UBLAS_INLINE
2860	bool same_closure (const matrix_binary_scalar2 &mbs2) const {
2861	return (*this).e1_.same_closure (mbs2.e1_) &&
2862	&e2_ == &(mbs2.e2_);
2863	}
2864
2865	// Iterator types
2866	private:
2867	typedef typename E1::const_iterator1 const_iterator11_type;
2868	typedef typename E1::const_iterator2 const_iterator12_type;
2869	typedef expression2_type const_subiterator2_type;
2870	typedef const value_type *const_pointer;
2871
2872	public:
2873	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2874	typedef indexed_const_iterator1<const_closure_type, typename const_iterator11_type::iterator_category> const_iterator1;
2875	typedef const_iterator1 iterator1;
2876	typedef indexed_const_iterator2<const_closure_type, typename const_iterator12_type::iterator_category> const_iterator2;
2877	typedef const_iterator2 iterator2;
2878	#else
2879	class const_iterator1;
2880	typedef const_iterator1 iterator1;
2881	class const_iterator2;
2882	typedef const_iterator2 iterator2;
2883	#endif
2884	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
2885	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
2886
2887	// Element lookup
2888	BOOST_UBLAS_INLINE
2889	const_iterator1 find1 (int rank, size_type i, size_type j) const {
2890	const_iterator11_type it11 (e1_.find1 (rank, i, j));
2891	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2892	return const_iterator1 (*this, it11.index1 (), it11.index2 ());
2893	#else
2894	return const_iterator1 (*this, it11, const_subiterator2_type (e2_));
2895	#endif
2896	}
2897	BOOST_UBLAS_INLINE
2898	const_iterator2 find2 (int rank, size_type i, size_type j) const {
2899	const_iterator12_type it12 (e1_.find2 (rank, i, j));
2900	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2901	return const_iterator2 (*this, it12.index1 (), it12.index2 ());
2902	#else
2903	return const_iterator2 (*this, it12, const_subiterator2_type (e2_));
2904	#endif
2905	}
2906
2907	// Iterators enhance the iterators of the referenced expression
2908	// with the binary functor.
2909
2910	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2911	class const_iterator1:
2912	public container_const_reference<matrix_binary_scalar2>,
2913	public iterator_base_traits<typename E1::const_iterator1::iterator_category>::template
2914	iterator_base<const_iterator1, value_type>::type {
2915	public:
2916	typedef typename E1::const_iterator1::iterator_category iterator_category;
2917	typedef typename matrix_binary_scalar2::difference_type difference_type;
2918	typedef typename matrix_binary_scalar2::value_type value_type;
2919	typedef typename matrix_binary_scalar2::const_reference reference;
2920	typedef typename matrix_binary_scalar2::const_pointer pointer;
2921
2922	typedef const_iterator2 dual_iterator_type;
2923	typedef const_reverse_iterator2 dual_reverse_iterator_type;
2924
2925	// Construction and destruction
2926	BOOST_UBLAS_INLINE
2927	const_iterator1 ():
2928	container_const_reference<self_type> (), it1_ (), it2_ () {}
2929	BOOST_UBLAS_INLINE
2930	const_iterator1 (const self_type &mbs, const const_iterator11_type &it1, const const_subiterator2_type &it2):
2931	container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
2932
2933	// Arithmetic
2934	BOOST_UBLAS_INLINE
2935	const_iterator1 &operator ++ () {
2936	++ it1_;
2937	return *this;
2938	}
2939	BOOST_UBLAS_INLINE
2940	const_iterator1 &operator -- () {
2941	-- it1_ ;
2942	return *this;
2943	}
2944	BOOST_UBLAS_INLINE
2945	const_iterator1 &operator += (difference_type n) {
2946	it1_ += n;
2947	return *this;
2948	}
2949	BOOST_UBLAS_INLINE
2950	const_iterator1 &operator -= (difference_type n) {
2951	it1_ -= n;
2952	return *this;
2953	}
2954	BOOST_UBLAS_INLINE
2955	difference_type operator - (const const_iterator1 &it) const {
2956	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2957	// FIXME we shouldn't compare floats
2958	// BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
2959	return it1_ - it.it1_;
2960	}
2961
2962	// Dereference
2963	BOOST_UBLAS_INLINE
2964	const_reference operator * () const {
2965	return functor_type::apply (*it1_, it2_);
2966	}
2967
2968	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2969	BOOST_UBLAS_INLINE
2970	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2971	typename self_type::
2972	#endif
2973	const_iterator2 begin () const {
2974	return (*this) ().find2 (1, index1 (), 0);
2975	}
2976	BOOST_UBLAS_INLINE
2977	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2978	typename self_type::
2979	#endif
2980	const_iterator2 end () const {
2981	return (this) ().find2 (1, index1 (), (this) ().size2 ());
2982	}
2983	BOOST_UBLAS_INLINE
2984	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2985	typename self_type::
2986	#endif
2987	const_reverse_iterator2 rbegin () const {
2988	return const_reverse_iterator2 (end ());
2989	}
2990	BOOST_UBLAS_INLINE
2991	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2992	typename self_type::
2993	#endif
2994	const_reverse_iterator2 rend () const {
2995	return const_reverse_iterator2 (begin ());
2996	}
2997	#endif
2998
2999	// Indices
3000	BOOST_UBLAS_INLINE
3001	size_type index1 () const {
3002	return it1_.index1 ();
3003	}
3004	BOOST_UBLAS_INLINE
3005	size_type index2 () const {
3006	return it1_.index2 ();
3007	}
3008
3009	// Assignment
3010	BOOST_UBLAS_INLINE
3011	const_iterator1 &operator = (const const_iterator1 &it) {
3012	container_const_reference<self_type>::assign (&it ());
3013	it1_ = it.it1_;
3014	it2_ = it.it2_;
3015	return *this;
3016	}
3017
3018	// Comparison
3019	BOOST_UBLAS_INLINE
3020	bool operator == (const const_iterator1 &it) const {
3021	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3022	// FIXME we shouldn't compare floats
3023	// BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3024	return it1_ == it.it1_;
3025	}
3026	BOOST_UBLAS_INLINE
3027	bool operator < (const const_iterator1 &it) const {
3028	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3029	// FIXME we shouldn't compare floats
3030	// BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3031	return it1_ < it.it1_;
3032	}
3033
3034	private:
3035	const_iterator11_type it1_;
3036	const_subiterator2_type it2_;
3037	};
3038	#endif
3039
3040	BOOST_UBLAS_INLINE
3041	const_iterator1 begin1 () const {
3042	return find1 (0, 0, 0);
3043	}
3044	BOOST_UBLAS_INLINE
3045	const_iterator1 end1 () const {
3046	return find1 (0, size1 (), 0);
3047	}
3048
3049	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3050	class const_iterator2:
3051	public container_const_reference<matrix_binary_scalar2>,
3052	public iterator_base_traits<typename E1::const_iterator2::iterator_category>::template
3053	iterator_base<const_iterator2, value_type>::type {
3054	public:
3055	typedef typename E1::const_iterator2::iterator_category iterator_category;
3056	typedef typename matrix_binary_scalar2::difference_type difference_type;
3057	typedef typename matrix_binary_scalar2::value_type value_type;
3058	typedef typename matrix_binary_scalar2::const_reference reference;
3059	typedef typename matrix_binary_scalar2::const_pointer pointer;
3060
3061	typedef const_iterator1 dual_iterator_type;
3062	typedef const_reverse_iterator1 dual_reverse_iterator_type;
3063
3064	// Construction and destruction
3065	BOOST_UBLAS_INLINE
3066	const_iterator2 ():
3067	container_const_reference<self_type> (), it1_ (), it2_ () {}
3068	BOOST_UBLAS_INLINE
3069	const_iterator2 (const self_type &mbs, const const_iterator12_type &it1, const const_subiterator2_type &it2):
3070	container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3071
3072	// Arithmetic
3073	BOOST_UBLAS_INLINE
3074	const_iterator2 &operator ++ () {
3075	++ it1_;
3076	return *this;
3077	}
3078	BOOST_UBLAS_INLINE
3079	const_iterator2 &operator -- () {
3080	-- it1_;
3081	return *this;
3082	}
3083	BOOST_UBLAS_INLINE
3084	const_iterator2 &operator += (difference_type n) {
3085	it1_ += n;
3086	return *this;
3087	}
3088	BOOST_UBLAS_INLINE
3089	const_iterator2 &operator -= (difference_type n) {
3090	it1_ -= n;
3091	return *this;
3092	}
3093	BOOST_UBLAS_INLINE
3094	difference_type operator - (const const_iterator2 &it) const {
3095	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3096	// FIXME we shouldn't compare floats
3097	// BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3098	return it1_ - it.it1_;
3099	}
3100
3101	// Dereference
3102	BOOST_UBLAS_INLINE
3103	const_reference operator * () const {
3104	return functor_type::apply (*it1_, it2_);
3105	}
3106
3107	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3108	BOOST_UBLAS_INLINE
3109	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3110	typename self_type::
3111	#endif
3112	const_iterator1 begin () const {
3113	return (*this) ().find1 (1, 0, index2 ());
3114	}
3115	BOOST_UBLAS_INLINE
3116	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3117	typename self_type::
3118	#endif
3119	const_iterator1 end () const {
3120	return (this) ().find1 (1, (this) ().size1 (), index2 ());
3121	}
3122	BOOST_UBLAS_INLINE
3123	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3124	typename self_type::
3125	#endif
3126	const_reverse_iterator1 rbegin () const {
3127	return const_reverse_iterator1 (end ());
3128	}
3129	BOOST_UBLAS_INLINE
3130	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3131	typename self_type::
3132	#endif
3133	const_reverse_iterator1 rend () const {
3134	return const_reverse_iterator1 (begin ());
3135	}
3136	#endif
3137
3138	// Indices
3139	BOOST_UBLAS_INLINE
3140	size_type index1 () const {
3141	return it1_.index1 ();
3142	}
3143	BOOST_UBLAS_INLINE
3144	size_type index2 () const {
3145	return it1_.index2 ();
3146	}
3147
3148	// Assignment
3149	BOOST_UBLAS_INLINE
3150	const_iterator2 &operator = (const const_iterator2 &it) {
3151	container_const_reference<self_type>::assign (&it ());
3152	it1_ = it.it1_;
3153	it2_ = it.it2_;
3154	return *this;
3155	}
3156
3157	// Comparison
3158	BOOST_UBLAS_INLINE
3159	bool operator == (const const_iterator2 &it) const {
3160	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3161	// FIXME we shouldn't compare floats
3162	// BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3163	return it1_ == it.it1_;
3164	}
3165	BOOST_UBLAS_INLINE
3166	bool operator < (const const_iterator2 &it) const {
3167	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3168	// FIXME we shouldn't compare floats
3169	// BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3170	return it1_ < it.it1_;
3171	}
3172
3173	private:
3174	const_iterator12_type it1_;
3175	const_subiterator2_type it2_;
3176	};
3177	#endif
3178
3179	BOOST_UBLAS_INLINE
3180	const_iterator2 begin2 () const {
3181	return find2 (0, 0, 0);
3182	}
3183	BOOST_UBLAS_INLINE
3184	const_iterator2 end2 () const {
3185	return find2 (0, 0, size2 ());
3186	}
3187
3188	// Reverse iterators
3189
3190	BOOST_UBLAS_INLINE
3191	const_reverse_iterator1 rbegin1 () const {
3192	return const_reverse_iterator1 (end1 ());
3193	}
3194	BOOST_UBLAS_INLINE
3195	const_reverse_iterator1 rend1 () const {
3196	return const_reverse_iterator1 (begin1 ());
3197	}
3198
3199	BOOST_UBLAS_INLINE
3200	const_reverse_iterator2 rbegin2 () const {
3201	return const_reverse_iterator2 (end2 ());
3202	}
3203	BOOST_UBLAS_INLINE
3204	const_reverse_iterator2 rend2 () const {
3205	return const_reverse_iterator2 (begin2 ());
3206	}
3207
3208	private:
3209	expression1_closure_type e1_;
3210	expression2_closure_type e2_;
3211	};
3212
3213	template<class E1, class E2, class F>
3214	struct matrix_binary_scalar2_traits {
3215	typedef matrix_binary_scalar2<E1, E2, F> expression_type; // allow E2 to be builtin type
3216	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
3217	typedef expression_type result_type;
3218	#else
3219	typedef typename E1::matrix_temporary_type result_type;
3220	#endif
3221	};
3222
3223	// (m * t) [i] [j] = m [i] [j] * t
3224	template<class E1, class T2>
3225	BOOST_UBLAS_INLINE
3226	typename matrix_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::result_type
3227	operator * (const matrix_expression<E1> &e1,
3228	const T2 &e2) {
3229	typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::expression_type expression_type;
3230	return expression_type (e1 (), e2);
3231	}
3232
3233	// (m / t) [i] [j] = m [i] [j] / t
3234	template<class E1, class T2>
3235	BOOST_UBLAS_INLINE
3236	typename matrix_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::result_type
3237	operator / (const matrix_expression<E1> &e1,
3238	const T2 &e2) {
3239	typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::expression_type expression_type;
3240	return expression_type (e1 (), e2);
3241	}
3242
3243
3244	template<class E1, class E2, class F>
3245	class matrix_vector_binary1:
3246	public vector_expression<matrix_vector_binary1<E1, E2, F> > {
3247
3248	public:
3249	typedef E1 expression1_type;
3250	typedef E2 expression2_type;
3251	private:
3252	typedef F functor_type;
3253	public:
3254	typedef typename E1::const_closure_type expression1_closure_type;
3255	typedef typename E2::const_closure_type expression2_closure_type;
3256	private:
3257	typedef matrix_vector_binary1<E1, E2, F> self_type;
3258	public:
3259	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
3260	using vector_expression<self_type>::operator ();
3261	#endif
3262	static const unsigned complexity = 1;
3263	typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
3264	typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
3265	typedef typename F::result_type value_type;
3266	typedef value_type const_reference;
3267	typedef const_reference reference;
3268	typedef const self_type const_closure_type;
3269	typedef const_closure_type closure_type;
3270	typedef unknown_storage_tag storage_category;
3271
3272	// Construction and destruction
3273	BOOST_UBLAS_INLINE
3274	matrix_vector_binary1 (const expression1_type &e1, const expression2_type &e2):
3275	e1_ (e1), e2_ (e2) {}
3276
3277	// Accessors
3278	BOOST_UBLAS_INLINE
3279	size_type size () const {
3280	return e1_.size1 ();
3281	}
3282
3283	public:
3284	// Expression accessors
3285	BOOST_UBLAS_INLINE
3286	const expression1_closure_type &expression1 () const {
3287	return e1_;
3288	}
3289	BOOST_UBLAS_INLINE
3290	const expression2_closure_type &expression2 () const {
3291	return e2_;
3292	}
3293
3294	public:
3295	// Element access
3296	BOOST_UBLAS_INLINE
3297	const_reference operator () (size_type i) const {
3298	return functor_type::apply (e1_, e2_, i);
3299	}
3300
3301	// Closure comparison
3302	BOOST_UBLAS_INLINE
3303	bool same_closure (const matrix_vector_binary1 &mvb1) const {
3304	return (*this).expression1 ().same_closure (mvb1.expression1 ()) &&
3305	(*this).expression2 ().same_closure (mvb1.expression2 ());
3306	}
3307
3308	// Iterator types
3309	private:
3310	typedef typename E1::const_iterator1 const_subiterator1_type;
3311	typedef typename E2::const_iterator const_subiterator2_type;
3312	typedef const value_type *const_pointer;
3313
3314	public:
3315	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3316	typedef indexed_const_iterator<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator;
3317	typedef const_iterator iterator;
3318	#else
3319	class const_iterator;
3320	typedef const_iterator iterator;
3321	#endif
3322
3323	// Element lookup
3324	BOOST_UBLAS_INLINE
3325	const_iterator find (size_type i) const {
3326	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3327	const_subiterator1_type it1 (e1_.find1 (0, i, 0));
3328	return const_iterator (*this, it1.index1 ());
3329	#else
3330	return const_iterator (*this, e1_.find1 (0, i, 0));
3331	#endif
3332	}
3333
3334
3335	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3336	class const_iterator:
3337	public container_const_reference<matrix_vector_binary1>,
3338	public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
3339	typename E2::const_iterator::iterator_category>::iterator_category>::template
3340	iterator_base<const_iterator, value_type>::type {
3341	public:
3342	typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
3343	typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
3344	typedef typename matrix_vector_binary1::difference_type difference_type;
3345	typedef typename matrix_vector_binary1::value_type value_type;
3346	typedef typename matrix_vector_binary1::const_reference reference;
3347	typedef typename matrix_vector_binary1::const_pointer pointer;
3348
3349	// Construction and destruction
3350	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3351	BOOST_UBLAS_INLINE
3352	const_iterator ():
3353	container_const_reference<self_type> (), it1_ (), e2_begin_ (), e2_end_ () {}
3354	BOOST_UBLAS_INLINE
3355	const_iterator (const self_type &mvb, const const_subiterator1_type &it1):
3356	container_const_reference<self_type> (mvb), it1_ (it1), e2_begin_ (mvb.expression2 ().begin ()), e2_end_ (mvb.expression2 ().end ()) {}
3357	#else
3358	BOOST_UBLAS_INLINE
3359	const_iterator ():
3360	container_const_reference<self_type> (), it1_ () {}
3361	BOOST_UBLAS_INLINE
3362	const_iterator (const self_type &mvb, const const_subiterator1_type &it1):
3363	container_const_reference<self_type> (mvb), it1_ (it1) {}
3364	#endif
3365
3366	// Dense random access specialization
3367	BOOST_UBLAS_INLINE
3368	value_type dereference (dense_random_access_iterator_tag) const {
3369	const self_type &mvb = (*this) ();
3370	#ifdef BOOST_UBLAS_USE_INDEXING
3371	return mvb (index ());
3372	#elif BOOST_UBLAS_USE_ITERATING
3373	difference_type size = BOOST_UBLAS_SAME (mvb.expression1 ().size2 (), mvb.expression2 ().size ());
3374	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3375	return functor_type::apply (size, it1_.begin (), e2_begin_);
3376	#else
3377	return functor_type::apply (size, it1_.begin (), mvb.expression2 ().begin ());
3378	#endif
3379	#else
3380	difference_type size = BOOST_UBLAS_SAME (mvb.expression1 ().size2 (), mvb.expression2 ().size ());
3381	if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
3382	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3383	return functor_type::apply (size, it1_.begin (), e2_begin_);
3384	#else
3385	return functor_type::apply (size, it1_.begin (), mvb.expression2 ().begin ());
3386	#endif
3387	else
3388	return mvb (index ());
3389	#endif
3390	}
3391
3392	// Packed bidirectional specialization
3393	BOOST_UBLAS_INLINE
3394	value_type dereference (packed_random_access_iterator_tag) const {
3395	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3396	return functor_type::apply (it1_.begin (), it1_.end (), e2_begin_, e2_end_);
3397	#else
3398	const self_type &mvb = (*this) ();
3399	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3400	return functor_type::apply (it1_.begin (), it1_.end (),
3401	mvb.expression2 ().begin (), mvb.expression2 ().end ());
3402	#else
3403	return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
3404	boost::numeric::ublas::end (it1_, iterator1_tag ()),
3405	mvb.expression2 ().begin (), mvb.expression2 ().end ());
3406	#endif
3407	#endif
3408	}
3409
3410	// Sparse bidirectional specialization
3411	BOOST_UBLAS_INLINE
3412	value_type dereference (sparse_bidirectional_iterator_tag) const {
3413	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3414	return functor_type::apply (it1_.begin (), it1_.end (), e2_begin_, e2_end_, sparse_bidirectional_iterator_tag ());
3415	#else
3416	const self_type &mvb = (*this) ();
3417	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3418	return functor_type::apply (it1_.begin (), it1_.end (),
3419	mvb.expression2 ().begin (), mvb.expression2 ().end (), sparse_bidirectional_iterator_tag ());
3420	#else
3421	return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
3422	boost::numeric::ublas::end (it1_, iterator1_tag ()),
3423	mvb.expression2 ().begin (), mvb.expression2 ().end (), sparse_bidirectional_iterator_tag ());
3424	#endif
3425	#endif
3426	}
3427
3428	// Arithmetic
3429	BOOST_UBLAS_INLINE
3430	const_iterator &operator ++ () {
3431	++ it1_;
3432	return *this;
3433	}
3434	BOOST_UBLAS_INLINE
3435	const_iterator &operator -- () {
3436	-- it1_;
3437	return *this;
3438	}
3439	BOOST_UBLAS_INLINE
3440	const_iterator &operator += (difference_type n) {
3441	it1_ += n;
3442	return *this;
3443	}
3444	BOOST_UBLAS_INLINE
3445	const_iterator &operator -= (difference_type n) {
3446	it1_ -= n;
3447	return *this;
3448	}
3449	BOOST_UBLAS_INLINE
3450	difference_type operator - (const const_iterator &it) const {
3451	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3452	return it1_ - it.it1_;
3453	}
3454
3455	// Dereference
3456	BOOST_UBLAS_INLINE
3457	const_reference operator * () const {
3458	return dereference (iterator_category ());
3459	}
3460
3461	// Index
3462	BOOST_UBLAS_INLINE
3463	size_type index () const {
3464	return it1_.index1 ();
3465	}
3466
3467	// Assignment
3468	BOOST_UBLAS_INLINE
3469	const_iterator &operator = (const const_iterator &it) {
3470	container_const_reference<self_type>::assign (&it ());
3471	it1_ = it.it1_;
3472	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3473	e2_begin_ = it.e2_begin_;
3474	e2_end_ = it.e2_end_;
3475	#endif
3476	return *this;
3477	}
3478
3479	// Comparison
3480	BOOST_UBLAS_INLINE
3481	bool operator == (const const_iterator &it) const {
3482	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3483	return it1_ == it.it1_;
3484	}
3485	BOOST_UBLAS_INLINE
3486	bool operator < (const const_iterator &it) const {
3487	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3488	return it1_ < it.it1_;
3489	}
3490
3491	private:
3492	const_subiterator1_type it1_;
3493	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3494	// Mutable due to assignment
3495	/* const */ const_subiterator2_type e2_begin_;
3496	/* const */ const_subiterator2_type e2_end_;
3497	#endif
3498	};
3499	#endif
3500
3501	BOOST_UBLAS_INLINE
3502	const_iterator begin () const {
3503	return find (0);
3504	}
3505	BOOST_UBLAS_INLINE
3506	const_iterator end () const {
3507	return find (size ());
3508	}
3509
3510	// Reverse iterator
3511	typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
3512
3513	BOOST_UBLAS_INLINE
3514	const_reverse_iterator rbegin () const {
3515	return const_reverse_iterator (end ());
3516	}
3517	BOOST_UBLAS_INLINE
3518	const_reverse_iterator rend () const {
3519	return const_reverse_iterator (begin ());
3520	}
3521
3522	private:
3523	expression1_closure_type e1_;
3524	expression2_closure_type e2_;
3525	};
3526
3527	template<class T1, class E1, class T2, class E2>
3528	struct matrix_vector_binary1_traits {
3529	typedef unknown_storage_tag storage_category;
3530	typedef row_major_tag orientation_category;
3531	typedef typename promote_traits<T1, T2>::promote_type promote_type;
3532	typedef matrix_vector_binary1<E1, E2, matrix_vector_prod1<T1, T2, promote_type> > expression_type;
3533	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
3534	typedef expression_type result_type;
3535	#else
3536	typedef typename E1::vector_temporary_type result_type;
3537	#endif
3538	};
3539
3540	template<class E1, class E2>
3541	BOOST_UBLAS_INLINE
3542	typename matrix_vector_binary1_traits<typename E1::value_type, E1,
3543	typename E2::value_type, E2>::result_type
3544	prod (const matrix_expression<E1> &e1,
3545	const vector_expression<E2> &e2,
3546	unknown_storage_tag,
3547	row_major_tag) {
3548	typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
3549	typename E2::value_type, E2>::expression_type expression_type;
3550	return expression_type (e1 (), e2 ());
3551	}
3552
3553	// Dispatcher
3554	template<class E1, class E2>
3555	BOOST_UBLAS_INLINE
3556	typename matrix_vector_binary1_traits<typename E1::value_type, E1,
3557	typename E2::value_type, E2>::result_type
3558	prod (const matrix_expression<E1> &e1,
3559	const vector_expression<E2> &e2) {
3560	BOOST_STATIC_ASSERT (E2::complexity == 0);
3561	typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
3562	typename E2::value_type, E2>::storage_category storage_category;
3563	typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
3564	typename E2::value_type, E2>::orientation_category orientation_category;
3565	return prod (e1, e2, storage_category (), orientation_category ());
3566	}
3567
3568	template<class E1, class E2>
3569	BOOST_UBLAS_INLINE
3570	typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3571	typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
3572	prec_prod (const matrix_expression<E1> &e1,
3573	const vector_expression<E2> &e2,
3574	unknown_storage_tag,
3575	row_major_tag) {
3576	typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3577	typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
3578	return expression_type (e1 (), e2 ());
3579	}
3580
3581	// Dispatcher
3582	template<class E1, class E2>
3583	BOOST_UBLAS_INLINE
3584	typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3585	typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
3586	prec_prod (const matrix_expression<E1> &e1,
3587	const vector_expression<E2> &e2) {
3588	BOOST_STATIC_ASSERT (E2::complexity == 0);
3589	typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3590	typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
3591	typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3592	typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
3593	return prec_prod (e1, e2, storage_category (), orientation_category ());
3594	}
3595
3596	template<class V, class E1, class E2>
3597	BOOST_UBLAS_INLINE
3598	V &
3599	prod (const matrix_expression<E1> &e1,
3600	const vector_expression<E2> &e2,
3601	V &v) {
3602	return v.assign (prod (e1, e2));
3603	}
3604
3605	template<class V, class E1, class E2>
3606	BOOST_UBLAS_INLINE
3607	V &
3608	prec_prod (const matrix_expression<E1> &e1,
3609	const vector_expression<E2> &e2,
3610	V &v) {
3611	return v.assign (prec_prod (e1, e2));
3612	}
3613
3614	template<class V, class E1, class E2>
3615	BOOST_UBLAS_INLINE
3616	V
3617	prod (const matrix_expression<E1> &e1,
3618	const vector_expression<E2> &e2) {
3619	return V (prod (e1, e2));
3620	}
3621
3622	template<class V, class E1, class E2>
3623	BOOST_UBLAS_INLINE
3624	V
3625	prec_prod (const matrix_expression<E1> &e1,
3626	const vector_expression<E2> &e2) {
3627	return V (prec_prod (e1, e2));
3628	}
3629
3630	template<class E1, class E2, class F>
3631	class matrix_vector_binary2:
3632	public vector_expression<matrix_vector_binary2<E1, E2, F> > {
3633
3634	typedef E1 expression1_type;
3635	typedef E2 expression2_type;
3636	typedef F functor_type;
3637	public:
3638	typedef typename E1::const_closure_type expression1_closure_type;
3639	typedef typename E2::const_closure_type expression2_closure_type;
3640	private:
3641	typedef matrix_vector_binary2<E1, E2, F> self_type;
3642	public:
3643	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
3644	using vector_expression<self_type>::operator ();
3645	#endif
3646	static const unsigned complexity = 1;
3647	typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
3648	typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
3649	typedef typename F::result_type value_type;
3650	typedef value_type const_reference;
3651	typedef const_reference reference;
3652	typedef const self_type const_closure_type;
3653	typedef const_closure_type closure_type;
3654	typedef unknown_storage_tag storage_category;
3655
3656	// Construction and destruction
3657	BOOST_UBLAS_INLINE
3658	matrix_vector_binary2 (const expression1_type &e1, const expression2_type &e2):
3659	e1_ (e1), e2_ (e2) {}
3660
3661	// Accessors
3662	BOOST_UBLAS_INLINE
3663	size_type size () const {
3664	return e2_.size2 ();
3665	}
3666
3667	public:
3668	// Expression accessors
3669	BOOST_UBLAS_INLINE
3670	const expression1_closure_type &expression1 () const {
3671	return e1_;
3672	}
3673	BOOST_UBLAS_INLINE
3674	const expression2_closure_type &expression2 () const {
3675	return e2_;
3676	}
3677	public:
3678
3679	// Element access
3680	BOOST_UBLAS_INLINE
3681	const_reference operator () (size_type j) const {
3682	return functor_type::apply (e1_, e2_, j);
3683	}
3684
3685	// Closure comparison
3686	BOOST_UBLAS_INLINE
3687	bool same_closure (const matrix_vector_binary2 &mvb2) const {
3688	return (*this).expression1 ().same_closure (mvb2.expression1 ()) &&
3689	(*this).expression2 ().same_closure (mvb2.expression2 ());
3690	}
3691
3692	// Iterator types
3693	private:
3694	typedef typename E1::const_iterator const_subiterator1_type;
3695	typedef typename E2::const_iterator2 const_subiterator2_type;
3696	typedef const value_type *const_pointer;
3697
3698	public:
3699	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3700	typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
3701	typedef const_iterator iterator;
3702	#else
3703	class const_iterator;
3704	typedef const_iterator iterator;
3705	#endif
3706
3707	// Element lookup
3708	BOOST_UBLAS_INLINE
3709	const_iterator find (size_type j) const {
3710	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3711	const_subiterator2_type it2 (e2_.find2 (0, 0, j));
3712	return const_iterator (*this, it2.index2 ());
3713	#else
3714	return const_iterator (*this, e2_.find2 (0, 0, j));
3715	#endif
3716	}
3717
3718
3719	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3720	class const_iterator:
3721	public container_const_reference<matrix_vector_binary2>,
3722	public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
3723	typename E2::const_iterator2::iterator_category>::iterator_category>::template
3724	iterator_base<const_iterator, value_type>::type {
3725	public:
3726	typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
3727	typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
3728	typedef typename matrix_vector_binary2::difference_type difference_type;
3729	typedef typename matrix_vector_binary2::value_type value_type;
3730	typedef typename matrix_vector_binary2::const_reference reference;
3731	typedef typename matrix_vector_binary2::const_pointer pointer;
3732
3733	// Construction and destruction
3734	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3735	BOOST_UBLAS_INLINE
3736	const_iterator ():
3737	container_const_reference<self_type> (), it2_ (), e1_begin_ (), e1_end_ () {}
3738	BOOST_UBLAS_INLINE
3739	const_iterator (const self_type &mvb, const const_subiterator2_type &it2):
3740	container_const_reference<self_type> (mvb), it2_ (it2), e1_begin_ (mvb.expression1 ().begin ()), e1_end_ (mvb.expression1 ().end ()) {}
3741	#else
3742	BOOST_UBLAS_INLINE
3743	const_iterator ():
3744	container_const_reference<self_type> (), it2_ () {}
3745	BOOST_UBLAS_INLINE
3746	const_iterator (const self_type &mvb, const const_subiterator2_type &it2):
3747	container_const_reference<self_type> (mvb), it2_ (it2) {}
3748	#endif
3749
3750	// Dense random access specialization
3751	BOOST_UBLAS_INLINE
3752	value_type dereference (dense_random_access_iterator_tag) const {
3753	const self_type &mvb = (*this) ();
3754	#ifdef BOOST_UBLAS_USE_INDEXING
3755	return mvb (index ());
3756	#elif BOOST_UBLAS_USE_ITERATING
3757	difference_type size = BOOST_UBLAS_SAME (mvb.expression2 ().size1 (), mvb.expression1 ().size ());
3758	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3759	return functor_type::apply (size, e1_begin_, it2_.begin ());
3760	#else
3761	return functor_type::apply (size, mvb.expression1 ().begin (), it2_.begin ());
3762	#endif
3763	#else
3764	difference_type size = BOOST_UBLAS_SAME (mvb.expression2 ().size1 (), mvb.expression1 ().size ());
3765	if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
3766	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3767	return functor_type::apply (size, e1_begin_, it2_.begin ());
3768	#else
3769	return functor_type::apply (size, mvb.expression1 ().begin (), it2_.begin ());
3770	#endif
3771	else
3772	return mvb (index ());
3773	#endif
3774	}
3775
3776	// Packed bidirectional specialization
3777	BOOST_UBLAS_INLINE
3778	value_type dereference (packed_random_access_iterator_tag) const {
3779	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3780	return functor_type::apply (e1_begin_, e1_end_, it2_.begin (), it2_.end ());
3781	#else
3782	const self_type &mvb = (*this) ();
3783	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3784	return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
3785	it2_.begin (), it2_.end ());
3786	#else
3787	return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
3788	boost::numeric::ublas::begin (it2_, iterator2_tag ()),
3789	boost::numeric::ublas::end (it2_, iterator2_tag ()));
3790	#endif
3791	#endif
3792	}
3793
3794	// Sparse bidirectional specialization
3795	BOOST_UBLAS_INLINE
3796	value_type dereference (sparse_bidirectional_iterator_tag) const {
3797	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3798	return functor_type::apply (e1_begin_, e1_end_, it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
3799	#else
3800	const self_type &mvb = (*this) ();
3801	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3802	return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
3803	it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
3804	#else
3805	return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
3806	boost::numeric::ublas::begin (it2_, iterator2_tag ()),
3807	boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
3808	#endif
3809	#endif
3810	}
3811
3812	// Arithmetic
3813	BOOST_UBLAS_INLINE
3814	const_iterator &operator ++ () {
3815	++ it2_;
3816	return *this;
3817	}
3818	BOOST_UBLAS_INLINE
3819	const_iterator &operator -- () {
3820	-- it2_;
3821	return *this;
3822	}
3823	BOOST_UBLAS_INLINE
3824	const_iterator &operator += (difference_type n) {
3825	it2_ += n;
3826	return *this;
3827	}
3828	BOOST_UBLAS_INLINE
3829	const_iterator &operator -= (difference_type n) {
3830	it2_ -= n;
3831	return *this;
3832	}
3833	BOOST_UBLAS_INLINE
3834	difference_type operator - (const const_iterator &it) const {
3835	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3836	return it2_ - it.it2_;
3837	}
3838
3839	// Dereference
3840	BOOST_UBLAS_INLINE
3841	const_reference operator * () const {
3842	return dereference (iterator_category ());
3843	}
3844
3845	// Index
3846	BOOST_UBLAS_INLINE
3847	size_type index () const {
3848	return it2_.index2 ();
3849	}
3850
3851	// Assignment
3852	BOOST_UBLAS_INLINE
3853	const_iterator &operator = (const const_iterator &it) {
3854	container_const_reference<self_type>::assign (&it ());
3855	it2_ = it.it2_;
3856	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3857	e1_begin_ = it.e1_begin_;
3858	e1_end_ = it.e1_end_;
3859	#endif
3860	return *this;
3861	}
3862
3863	// Comparison
3864	BOOST_UBLAS_INLINE
3865	bool operator == (const const_iterator &it) const {
3866	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3867	return it2_ == it.it2_;
3868	}
3869	BOOST_UBLAS_INLINE
3870	bool operator < (const const_iterator &it) const {
3871	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3872	return it2_ < it.it2_;
3873	}
3874
3875	private:
3876	const_subiterator2_type it2_;
3877	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
3878	// Mutable due to assignment
3879	/* const */ const_subiterator1_type e1_begin_;
3880	/* const */ const_subiterator1_type e1_end_;
3881	#endif
3882	};
3883	#endif
3884
3885	BOOST_UBLAS_INLINE
3886	const_iterator begin () const {
3887	return find (0);
3888	}
3889	BOOST_UBLAS_INLINE
3890	const_iterator end () const {
3891	return find (size ());
3892	}
3893
3894	// Reverse iterator
3895	typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
3896
3897	BOOST_UBLAS_INLINE
3898	const_reverse_iterator rbegin () const {
3899	return const_reverse_iterator (end ());
3900	}
3901	BOOST_UBLAS_INLINE
3902	const_reverse_iterator rend () const {
3903	return const_reverse_iterator (begin ());
3904	}
3905
3906	private:
3907	expression1_closure_type e1_;
3908	expression2_closure_type e2_;
3909	};
3910
3911	template<class T1, class E1, class T2, class E2>
3912	struct matrix_vector_binary2_traits {
3913	typedef unknown_storage_tag storage_category;
3914	typedef column_major_tag orientation_category;
3915	typedef typename promote_traits<T1, T2>::promote_type promote_type;
3916	typedef matrix_vector_binary2<E1, E2, matrix_vector_prod2<T1, T2, promote_type> > expression_type;
3917	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
3918	typedef expression_type result_type;
3919	#else
3920	typedef typename E2::vector_temporary_type result_type;
3921	#endif
3922	};
3923
3924	template<class E1, class E2>
3925	BOOST_UBLAS_INLINE
3926	typename matrix_vector_binary2_traits<typename E1::value_type, E1,
3927	typename E2::value_type, E2>::result_type
3928	prod (const vector_expression<E1> &e1,
3929	const matrix_expression<E2> &e2,
3930	unknown_storage_tag,
3931	column_major_tag) {
3932	typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
3933	typename E2::value_type, E2>::expression_type expression_type;
3934	return expression_type (e1 (), e2 ());
3935	}
3936
3937	// Dispatcher
3938	template<class E1, class E2>
3939	BOOST_UBLAS_INLINE
3940	typename matrix_vector_binary2_traits<typename E1::value_type, E1,
3941	typename E2::value_type, E2>::result_type
3942	prod (const vector_expression<E1> &e1,
3943	const matrix_expression<E2> &e2) {
3944	BOOST_STATIC_ASSERT (E1::complexity == 0);
3945	typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
3946	typename E2::value_type, E2>::storage_category storage_category;
3947	typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
3948	typename E2::value_type, E2>::orientation_category orientation_category;
3949	return prod (e1, e2, storage_category (), orientation_category ());
3950	}
3951
3952	template<class E1, class E2>
3953	BOOST_UBLAS_INLINE
3954	typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3955	typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
3956	prec_prod (const vector_expression<E1> &e1,
3957	const matrix_expression<E2> &e2,
3958	unknown_storage_tag,
3959	column_major_tag) {
3960	typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3961	typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
3962	return expression_type (e1 (), e2 ());
3963	}
3964
3965	// Dispatcher
3966	template<class E1, class E2>
3967	BOOST_UBLAS_INLINE
3968	typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3969	typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
3970	prec_prod (const vector_expression<E1> &e1,
3971	const matrix_expression<E2> &e2) {
3972	BOOST_STATIC_ASSERT (E1::complexity == 0);
3973	typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3974	typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
3975	typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
3976	typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
3977	return prec_prod (e1, e2, storage_category (), orientation_category ());
3978	}
3979
3980	template<class V, class E1, class E2>
3981	BOOST_UBLAS_INLINE
3982	V &
3983	prod (const vector_expression<E1> &e1,
3984	const matrix_expression<E2> &e2,
3985	V &v) {
3986	return v.assign (prod (e1, e2));
3987	}
3988
3989	template<class V, class E1, class E2>
3990	BOOST_UBLAS_INLINE
3991	V &
3992	prec_prod (const vector_expression<E1> &e1,
3993	const matrix_expression<E2> &e2,
3994	V &v) {
3995	return v.assign (prec_prod (e1, e2));
3996	}
3997
3998	template<class V, class E1, class E2>
3999	BOOST_UBLAS_INLINE
4000	V
4001	prod (const vector_expression<E1> &e1,
4002	const matrix_expression<E2> &e2) {
4003	return V (prod (e1, e2));
4004	}
4005
4006	template<class V, class E1, class E2>
4007	BOOST_UBLAS_INLINE
4008	V
4009	prec_prod (const vector_expression<E1> &e1,
4010	const matrix_expression<E2> &e2) {
4011	return V (prec_prod (e1, e2));
4012	}
4013
4014	template<class E1, class E2, class F>
4015	class matrix_matrix_binary:
4016	public matrix_expression<matrix_matrix_binary<E1, E2, F> > {
4017
4018	public:
4019	typedef E1 expression1_type;
4020	typedef E2 expression2_type;
4021	private:
4022	typedef F functor_type;
4023	public:
4024	typedef typename E1::const_closure_type expression1_closure_type;
4025	typedef typename E2::const_closure_type expression2_closure_type;
4026	private:
4027	typedef matrix_matrix_binary<E1, E2, F> self_type;
4028	public:
4029	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
4030	using matrix_expression<self_type>::operator ();
4031	#endif
4032	static const unsigned complexity = 1;
4033	typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
4034	typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
4035	typedef typename F::result_type value_type;
4036	typedef value_type const_reference;
4037	typedef const_reference reference;
4038	typedef const self_type const_closure_type;
4039	typedef const_closure_type closure_type;
4040	typedef unknown_orientation_tag orientation_category;
4041	typedef unknown_storage_tag storage_category;
4042
4043	// Construction and destruction
4044	BOOST_UBLAS_INLINE
4045	matrix_matrix_binary (const expression1_type &e1, const expression2_type &e2):
4046	e1_ (e1), e2_ (e2) {}
4047
4048	// Accessors
4049	BOOST_UBLAS_INLINE
4050	size_type size1 () const {
4051	return e1_.size1 ();
4052	}
4053	BOOST_UBLAS_INLINE
4054	size_type size2 () const {
4055	return e2_.size2 ();
4056	}
4057
4058	public:
4059	// Expression accessors
4060	BOOST_UBLAS_INLINE
4061	const expression1_closure_type &expression1 () const {
4062	return e1_;
4063	}
4064	BOOST_UBLAS_INLINE
4065	const expression2_closure_type &expression2 () const {
4066	return e2_;
4067	}
4068
4069	public:
4070	// Element access
4071	BOOST_UBLAS_INLINE
4072	const_reference operator () (size_type i, size_type j) const {
4073	return functor_type::apply (e1_, e2_, i, j);
4074	}
4075
4076	// Closure comparison
4077	BOOST_UBLAS_INLINE
4078	bool same_closure (const matrix_matrix_binary &mmb) const {
4079	return (*this).expression1 ().same_closure (mmb.expression1 ()) &&
4080	(*this).expression2 ().same_closure (mmb.expression2 ());
4081	}
4082
4083	// Iterator types
4084	private:
4085	typedef typename E1::const_iterator1 const_iterator11_type;
4086	typedef typename E1::const_iterator2 const_iterator12_type;
4087	typedef typename E2::const_iterator1 const_iterator21_type;
4088	typedef typename E2::const_iterator2 const_iterator22_type;
4089	typedef const value_type *const_pointer;
4090
4091	public:
4092	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4093	typedef typename iterator_restrict_traits<typename const_iterator11_type::iterator_category,
4094	typename const_iterator22_type::iterator_category>::iterator_category iterator_category;
4095	typedef indexed_const_iterator1<const_closure_type, iterator_category> const_iterator1;
4096	typedef const_iterator1 iterator1;
4097	typedef indexed_const_iterator2<const_closure_type, iterator_category> const_iterator2;
4098	typedef const_iterator2 iterator2;
4099	#else
4100	class const_iterator1;
4101	typedef const_iterator1 iterator1;
4102	class const_iterator2;
4103	typedef const_iterator2 iterator2;
4104	#endif
4105	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
4106	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
4107
4108	// Element lookup
4109	BOOST_UBLAS_INLINE
4110	const_iterator1 find1 (int /* rank */, size_type i, size_type j) const {
4111	// FIXME sparse matrix tests fail!
4112	// const_iterator11_type it11 (e1_.find1 (rank, i, 0));
4113	const_iterator11_type it11 (e1_.find1 (0, i, 0));
4114	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4115	return const_iterator1 (*this, it11.index1 (), j);
4116	#else
4117	// FIXME sparse matrix tests fail!
4118	// const_iterator22_type it22 (e2_.find2 (rank, 0, j));
4119	const_iterator22_type it22 (e2_.find2 (0, 0, j));
4120	return const_iterator1 (*this, it11, it22);
4121	#endif
4122	}
4123	BOOST_UBLAS_INLINE
4124	const_iterator2 find2 (int /* rank */, size_type i, size_type j) const {
4125	// FIXME sparse matrix tests fail!
4126	// const_iterator22_type it22 (e2_.find2 (rank, 0, j));
4127	const_iterator22_type it22 (e2_.find2 (0, 0, j));
4128	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4129	return const_iterator2 (*this, i, it22.index2 ());
4130	#else
4131	// FIXME sparse matrix tests fail!
4132	// const_iterator11_type it11 (e1_.find1 (rank, i, 0));
4133	const_iterator11_type it11 (e1_.find1 (0, i, 0));
4134	return const_iterator2 (*this, it11, it22);
4135	#endif
4136	}
4137
4138
4139	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
4140	class const_iterator1:
4141	public container_const_reference<matrix_matrix_binary>,
4142	public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4143	typename E2::const_iterator2::iterator_category>::iterator_category>::template
4144	iterator_base<const_iterator1, value_type>::type {
4145	public:
4146	typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4147	typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
4148	typedef typename matrix_matrix_binary::difference_type difference_type;
4149	typedef typename matrix_matrix_binary::value_type value_type;
4150	typedef typename matrix_matrix_binary::const_reference reference;
4151	typedef typename matrix_matrix_binary::const_pointer pointer;
4152
4153	typedef const_iterator2 dual_iterator_type;
4154	typedef const_reverse_iterator2 dual_reverse_iterator_type;
4155
4156	// Construction and destruction
4157	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4158	BOOST_UBLAS_INLINE
4159	const_iterator1 ():
4160	container_const_reference<self_type> (), it1_ (), it2_ (), it2_begin_ (), it2_end_ () {}
4161	BOOST_UBLAS_INLINE
4162	const_iterator1 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
4163	container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2), it2_begin_ (it2.begin ()), it2_end_ (it2.end ()) {}
4164	#else
4165	BOOST_UBLAS_INLINE
4166	const_iterator1 ():
4167	container_const_reference<self_type> (), it1_ (), it2_ () {}
4168	BOOST_UBLAS_INLINE
4169	const_iterator1 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
4170	container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2) {}
4171	#endif
4172
4173	// Random access specialization
4174	BOOST_UBLAS_INLINE
4175	value_type dereference (dense_random_access_iterator_tag) const {
4176	const self_type &mmb = (*this) ();
4177	#ifdef BOOST_UBLAS_USE_INDEXING
4178	return mmb (index1 (), index2 ());
4179	#elif BOOST_UBLAS_USE_ITERATING
4180	difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
4181	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4182	return functor_type::apply (size, it1_.begin (), it2_begin_);
4183	#else
4184	return functor_type::apply (size, it1_.begin (), it2_.begin ());
4185	#endif
4186	#else
4187	difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
4188	if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
4189	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4190	return functor_type::apply (size, it1_.begin (), it2_begin_);
4191	#else
4192	return functor_type::apply (size, it1_.begin (), it2_.begin ());
4193	#endif
4194	else
4195	return mmb (index1 (), index2 ());
4196	#endif
4197	}
4198
4199	// Packed bidirectional specialization
4200	BOOST_UBLAS_INLINE
4201	value_type dereference (packed_random_access_iterator_tag) const {
4202	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4203	return functor_type::apply (it1_.begin (), it1_.end (),
4204	it2_begin_, it2_end_, packed_random_access_iterator_tag ());
4205	#else
4206	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4207	return functor_type::apply (it1_.begin (), it1_.end (),
4208	it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
4209	#else
4210	return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4211	boost::numeric::ublas::end (it1_, iterator1_tag ()),
4212	boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4213	boost::numeric::ublas::end (it2_, iterator2_tag ()), packed_random_access_iterator_tag ());
4214	#endif
4215	#endif
4216	}
4217
4218	// Sparse bidirectional specialization
4219	BOOST_UBLAS_INLINE
4220	value_type dereference (sparse_bidirectional_iterator_tag) const {
4221	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4222	return functor_type::apply (it1_.begin (), it1_.end (),
4223	it2_begin_, it2_end_, sparse_bidirectional_iterator_tag ());
4224	#else
4225	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4226	return functor_type::apply (it1_.begin (), it1_.end (),
4227	it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
4228	#else
4229	return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4230	boost::numeric::ublas::end (it1_, iterator1_tag ()),
4231	boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4232	boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
4233	#endif
4234	#endif
4235	}
4236
4237	// Arithmetic
4238	BOOST_UBLAS_INLINE
4239	const_iterator1 &operator ++ () {
4240	++ it1_;
4241	return *this;
4242	}
4243	BOOST_UBLAS_INLINE
4244	const_iterator1 &operator -- () {
4245	-- it1_;
4246	return *this;
4247	}
4248	BOOST_UBLAS_INLINE
4249	const_iterator1 &operator += (difference_type n) {
4250	it1_ += n;
4251	return *this;
4252	}
4253	BOOST_UBLAS_INLINE
4254	const_iterator1 &operator -= (difference_type n) {
4255	it1_ -= n;
4256	return *this;
4257	}
4258	BOOST_UBLAS_INLINE
4259	difference_type operator - (const const_iterator1 &it) const {
4260	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4261	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
4262	return it1_ - it.it1_;
4263	}
4264
4265	// Dereference
4266	BOOST_UBLAS_INLINE
4267	const_reference operator * () const {
4268	return dereference (iterator_category ());
4269	}
4270
4271	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4272	BOOST_UBLAS_INLINE
4273	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4274	typename self_type::
4275	#endif
4276	const_iterator2 begin () const {
4277	return (*this) ().find2 (1, index1 (), 0);
4278	}
4279	BOOST_UBLAS_INLINE
4280	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4281	typename self_type::
4282	#endif
4283	const_iterator2 end () const {
4284	return (this) ().find2 (1, index1 (), (this) ().size2 ());
4285	}
4286	BOOST_UBLAS_INLINE
4287	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4288	typename self_type::
4289	#endif
4290	const_reverse_iterator2 rbegin () const {
4291	return const_reverse_iterator2 (end ());
4292	}
4293	BOOST_UBLAS_INLINE
4294	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4295	typename self_type::
4296	#endif
4297	const_reverse_iterator2 rend () const {
4298	return const_reverse_iterator2 (begin ());
4299	}
4300	#endif
4301
4302	// Indices
4303	BOOST_UBLAS_INLINE
4304	size_type index1 () const {
4305	return it1_.index1 ();
4306	}
4307	BOOST_UBLAS_INLINE
4308	size_type index2 () const {
4309	return it2_.index2 ();
4310	}
4311
4312	// Assignment
4313	BOOST_UBLAS_INLINE
4314	const_iterator1 &operator = (const const_iterator1 &it) {
4315	container_const_reference<self_type>::assign (&it ());
4316	it1_ = it.it1_;
4317	it2_ = it.it2_;
4318	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4319	it2_begin_ = it.it2_begin_;
4320	it2_end_ = it.it2_end_;
4321	#endif
4322	return *this;
4323	}
4324
4325	// Comparison
4326	BOOST_UBLAS_INLINE
4327	bool operator == (const const_iterator1 &it) const {
4328	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4329	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
4330	return it1_ == it.it1_;
4331	}
4332	BOOST_UBLAS_INLINE
4333	bool operator < (const const_iterator1 &it) const {
4334	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4335	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
4336	return it1_ < it.it1_;
4337	}
4338
4339	private:
4340	const_iterator11_type it1_;
4341	// Mutable due to assignment
4342	/* const */ const_iterator22_type it2_;
4343	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4344	/* const */ const_iterator21_type it2_begin_;
4345	/* const */ const_iterator21_type it2_end_;
4346	#endif
4347	};
4348	#endif
4349
4350	BOOST_UBLAS_INLINE
4351	const_iterator1 begin1 () const {
4352	return find1 (0, 0, 0);
4353	}
4354	BOOST_UBLAS_INLINE
4355	const_iterator1 end1 () const {
4356	return find1 (0, size1 (), 0);
4357	}
4358
4359	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
4360	class const_iterator2:
4361	public container_const_reference<matrix_matrix_binary>,
4362	public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4363	typename E2::const_iterator2::iterator_category>::iterator_category>::template
4364	iterator_base<const_iterator2, value_type>::type {
4365	public:
4366	typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4367	typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
4368	typedef typename matrix_matrix_binary::difference_type difference_type;
4369	typedef typename matrix_matrix_binary::value_type value_type;
4370	typedef typename matrix_matrix_binary::const_reference reference;
4371	typedef typename matrix_matrix_binary::const_pointer pointer;
4372
4373	typedef const_iterator1 dual_iterator_type;
4374	typedef const_reverse_iterator1 dual_reverse_iterator_type;
4375
4376	// Construction and destruction
4377	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4378	BOOST_UBLAS_INLINE
4379	const_iterator2 ():
4380	container_const_reference<self_type> (), it1_ (), it2_ (), it1_begin_ (), it1_end_ () {}
4381	BOOST_UBLAS_INLINE
4382	const_iterator2 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
4383	container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2), it1_begin_ (it1.begin ()), it1_end_ (it1.end ()) {}
4384	#else
4385	BOOST_UBLAS_INLINE
4386	const_iterator2 ():
4387	container_const_reference<self_type> (), it1_ (), it2_ () {}
4388	BOOST_UBLAS_INLINE
4389	const_iterator2 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
4390	container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2) {}
4391	#endif
4392
4393	// Random access specialization
4394	BOOST_UBLAS_INLINE
4395	value_type dereference (dense_random_access_iterator_tag) const {
4396	const self_type &mmb = (*this) ();
4397	#ifdef BOOST_UBLAS_USE_INDEXING
4398	return mmb (index1 (), index2 ());
4399	#elif BOOST_UBLAS_USE_ITERATING
4400	difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
4401	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4402	return functor_type::apply (size, it1_begin_, it2_.begin ());
4403	#else
4404	return functor_type::apply (size, it1_.begin (), it2_.begin ());
4405	#endif
4406	#else
4407	difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
4408	if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
4409	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4410	return functor_type::apply (size, it1_begin_, it2_.begin ());
4411	#else
4412	return functor_type::apply (size, it1_.begin (), it2_.begin ());
4413	#endif
4414	else
4415	return mmb (index1 (), index2 ());
4416	#endif
4417	}
4418
4419	// Packed bidirectional specialization
4420	BOOST_UBLAS_INLINE
4421	value_type dereference (packed_random_access_iterator_tag) const {
4422	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4423	return functor_type::apply (it1_begin_, it1_end_,
4424	it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
4425	#else
4426	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4427	return functor_type::apply (it1_.begin (), it1_.end (),
4428	it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
4429	#else
4430	return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4431	boost::numeric::ublas::end (it1_, iterator1_tag ()),
4432	boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4433	boost::numeric::ublas::end (it2_, iterator2_tag ()), packed_random_access_iterator_tag ());
4434	#endif
4435	#endif
4436	}
4437
4438	// Sparse bidirectional specialization
4439	BOOST_UBLAS_INLINE
4440	value_type dereference (sparse_bidirectional_iterator_tag) const {
4441	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4442	return functor_type::apply (it1_begin_, it1_end_,
4443	it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
4444	#else
4445	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4446	return functor_type::apply (it1_.begin (), it1_.end (),
4447	it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
4448	#else
4449	return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4450	boost::numeric::ublas::end (it1_, iterator1_tag ()),
4451	boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4452	boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
4453	#endif
4454	#endif
4455	}
4456
4457	// Arithmetic
4458	BOOST_UBLAS_INLINE
4459	const_iterator2 &operator ++ () {
4460	++ it2_;
4461	return *this;
4462	}
4463	BOOST_UBLAS_INLINE
4464	const_iterator2 &operator -- () {
4465	-- it2_;
4466	return *this;
4467	}
4468	BOOST_UBLAS_INLINE
4469	const_iterator2 &operator += (difference_type n) {
4470	it2_ += n;
4471	return *this;
4472	}
4473	BOOST_UBLAS_INLINE
4474	const_iterator2 &operator -= (difference_type n) {
4475	it2_ -= n;
4476	return *this;
4477	}
4478	BOOST_UBLAS_INLINE
4479	difference_type operator - (const const_iterator2 &it) const {
4480	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4481	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
4482	return it2_ - it.it2_;
4483	}
4484
4485	// Dereference
4486	BOOST_UBLAS_INLINE
4487	const_reference operator * () const {
4488	return dereference (iterator_category ());
4489	}
4490
4491	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4492	BOOST_UBLAS_INLINE
4493	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4494	typename self_type::
4495	#endif
4496	const_iterator1 begin () const {
4497	return (*this) ().find1 (1, 0, index2 ());
4498	}
4499	BOOST_UBLAS_INLINE
4500	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4501	typename self_type::
4502	#endif
4503	const_iterator1 end () const {
4504	return (this) ().find1 (1, (this) ().size1 (), index2 ());
4505	}
4506	BOOST_UBLAS_INLINE
4507	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4508	typename self_type::
4509	#endif
4510	const_reverse_iterator1 rbegin () const {
4511	return const_reverse_iterator1 (end ());
4512	}
4513	BOOST_UBLAS_INLINE
4514	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4515	typename self_type::
4516	#endif
4517	const_reverse_iterator1 rend () const {
4518	return const_reverse_iterator1 (begin ());
4519	}
4520	#endif
4521
4522	// Indices
4523	BOOST_UBLAS_INLINE
4524	size_type index1 () const {
4525	return it1_.index1 ();
4526	}
4527	BOOST_UBLAS_INLINE
4528	size_type index2 () const {
4529	return it2_.index2 ();
4530	}
4531
4532	// Assignment
4533	BOOST_UBLAS_INLINE
4534	const_iterator2 &operator = (const const_iterator2 &it) {
4535	container_const_reference<self_type>::assign (&it ());
4536	it1_ = it.it1_;
4537	it2_ = it.it2_;
4538	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4539	it1_begin_ = it.it1_begin_;
4540	it1_end_ = it.it1_end_;
4541	#endif
4542	return *this;
4543	}
4544
4545	// Comparison
4546	BOOST_UBLAS_INLINE
4547	bool operator == (const const_iterator2 &it) const {
4548	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4549	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
4550	return it2_ == it.it2_;
4551	}
4552	BOOST_UBLAS_INLINE
4553	bool operator < (const const_iterator2 &it) const {
4554	BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4555	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
4556	return it2_ < it.it2_;
4557	}
4558
4559	private:
4560	// Mutable due to assignment
4561	/* const */ const_iterator11_type it1_;
4562	const_iterator22_type it2_;
4563	#ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4564	/* const */ const_iterator12_type it1_begin_;
4565	/* const */ const_iterator12_type it1_end_;
4566	#endif
4567	};
4568	#endif
4569
4570	BOOST_UBLAS_INLINE
4571	const_iterator2 begin2 () const {
4572	return find2 (0, 0, 0);
4573	}
4574	BOOST_UBLAS_INLINE
4575	const_iterator2 end2 () const {
4576	return find2 (0, 0, size2 ());
4577	}
4578
4579	// Reverse iterators
4580
4581	BOOST_UBLAS_INLINE
4582	const_reverse_iterator1 rbegin1 () const {
4583	return const_reverse_iterator1 (end1 ());
4584	}
4585	BOOST_UBLAS_INLINE
4586	const_reverse_iterator1 rend1 () const {
4587	return const_reverse_iterator1 (begin1 ());
4588	}
4589
4590	BOOST_UBLAS_INLINE
4591	const_reverse_iterator2 rbegin2 () const {
4592	return const_reverse_iterator2 (end2 ());
4593	}
4594	BOOST_UBLAS_INLINE
4595	const_reverse_iterator2 rend2 () const {
4596	return const_reverse_iterator2 (begin2 ());
4597	}
4598
4599	private:
4600	expression1_closure_type e1_;
4601	expression2_closure_type e2_;
4602	};
4603
4604	template<class T1, class E1, class T2, class E2>
4605	struct matrix_matrix_binary_traits {
4606	typedef unknown_storage_tag storage_category;
4607	typedef unknown_orientation_tag orientation_category;
4608	typedef typename promote_traits<T1, T2>::promote_type promote_type;
4609	typedef matrix_matrix_binary<E1, E2, matrix_matrix_prod<T1, T2, promote_type> > expression_type;
4610	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
4611	typedef expression_type result_type;
4612	#else
4613	typedef typename E1::matrix_temporary_type result_type;
4614	#endif
4615	};
4616
4617	template<class E1, class E2>
4618	BOOST_UBLAS_INLINE
4619	typename matrix_matrix_binary_traits<typename E1::value_type, E1,
4620	typename E2::value_type, E2>::result_type
4621	prod (const matrix_expression<E1> &e1,
4622	const matrix_expression<E2> &e2,
4623	unknown_storage_tag,
4624	unknown_orientation_tag) {
4625	typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
4626	typename E2::value_type, E2>::expression_type expression_type;
4627	return expression_type (e1 (), e2 ());
4628	}
4629
4630	// Dispatcher
4631	template<class E1, class E2>
4632	BOOST_UBLAS_INLINE
4633	typename matrix_matrix_binary_traits<typename E1::value_type, E1,
4634	typename E2::value_type, E2>::result_type
4635	prod (const matrix_expression<E1> &e1,
4636	const matrix_expression<E2> &e2) {
4637	BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
4638	typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
4639	typename E2::value_type, E2>::storage_category storage_category;
4640	typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
4641	typename E2::value_type, E2>::orientation_category orientation_category;
4642	return prod (e1, e2, storage_category (), orientation_category ());
4643	}
4644
4645	template<class E1, class E2>
4646	BOOST_UBLAS_INLINE
4647	typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4648	typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4649	prec_prod (const matrix_expression<E1> &e1,
4650	const matrix_expression<E2> &e2,
4651	unknown_storage_tag,
4652	unknown_orientation_tag) {
4653	typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4654	typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
4655	return expression_type (e1 (), e2 ());
4656	}
4657
4658	// Dispatcher
4659	template<class E1, class E2>
4660	BOOST_UBLAS_INLINE
4661	typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4662	typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4663	prec_prod (const matrix_expression<E1> &e1,
4664	const matrix_expression<E2> &e2) {
4665	BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
4666	typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4667	typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
4668	typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4669	typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
4670	return prec_prod (e1, e2, storage_category (), orientation_category ());
4671	}
4672
4673	template<class M, class E1, class E2>
4674	BOOST_UBLAS_INLINE
4675	M &
4676	prod (const matrix_expression<E1> &e1,
4677	const matrix_expression<E2> &e2,
4678	M &m) {
4679	return m.assign (prod (e1, e2));
4680	}
4681
4682	template<class M, class E1, class E2>
4683	BOOST_UBLAS_INLINE
4684	M &
4685	prec_prod (const matrix_expression<E1> &e1,
4686	const matrix_expression<E2> &e2,
4687	M &m) {
4688	return m.assign (prec_prod (e1, e2));
4689	}
4690
4691	template<class M, class E1, class E2>
4692	BOOST_UBLAS_INLINE
4693	M
4694	prod (const matrix_expression<E1> &e1,
4695	const matrix_expression<E2> &e2) {
4696	return M (prod (e1, e2));
4697	}
4698
4699	template<class M, class E1, class E2>
4700	BOOST_UBLAS_INLINE
4701	M
4702	prec_prod (const matrix_expression<E1> &e1,
4703	const matrix_expression<E2> &e2) {
4704	return M (prec_prod (e1, e2));
4705	}
4706
4707	template<class E, class F>
4708	class matrix_scalar_unary:
4709	public scalar_expression<matrix_scalar_unary<E, F> > {
4710	public:
4711	typedef E expression_type;
4712	typedef F functor_type;
4713	typedef typename F::result_type value_type;
4714	typedef typename E::const_closure_type expression_closure_type;
4715
4716	// Construction and destruction
4717	BOOST_UBLAS_INLINE
4718	explicit matrix_scalar_unary (const expression_type &e):
4719	e_ (e) {}
4720
4721	private:
4722	// Expression accessors
4723	BOOST_UBLAS_INLINE
4724	const expression_closure_type &expression () const {
4725	return e_;
4726	}
4727
4728	public:
4729	BOOST_UBLAS_INLINE
4730	operator value_type () const {
4731	return functor_type::apply (e_);
4732	}
4733
4734	private:
4735	expression_closure_type e_;
4736	};
4737
4738	template<class E, class F>
4739	struct matrix_scalar_unary_traits {
4740	typedef matrix_scalar_unary<E, F> expression_type;
4741	#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
4742	typedef expression_type result_type;
4743	#else
4744	typedef typename F::result_type result_type;
4745	#endif
4746	};
4747
4748	template<class E>
4749	BOOST_UBLAS_INLINE
4750	typename matrix_scalar_unary_traits<E, matrix_norm_1<typename E::value_type> >::result_type
4751	norm_1 (const matrix_expression<E> &e) {
4752	typedef typename matrix_scalar_unary_traits<E, matrix_norm_1<typename E::value_type> >::expression_type expression_type;
4753	return expression_type (e ());
4754	}
4755
4756	template<class E>
4757	BOOST_UBLAS_INLINE
4758	typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<typename E::value_type> >::result_type
4759	norm_frobenius (const matrix_expression<E> &e) {
4760	typedef typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<typename E::value_type> >::expression_type expression_type;
4761	return expression_type (e ());
4762	}
4763
4764	template<class E>
4765	BOOST_UBLAS_INLINE
4766	typename matrix_scalar_unary_traits<E, matrix_norm_inf<typename E::value_type> >::result_type
4767	norm_inf (const matrix_expression<E> &e) {
4768	typedef typename matrix_scalar_unary_traits<E, matrix_norm_inf<typename E::value_type> >::expression_type expression_type;
4769	return expression_type (e ());
4770	}
4771
4772	}}}
4773
4774	#endif

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: