Planet

navi

home

PPS

about

screenshots

download

development

forum

Context Navigation

source: downloads/boost_1_34_1/boost/numeric/ublas/matrix_expression.hpp @ 33

Last change on this file since 33 was 29, checked in by landauf, 17 years ago
updated boost from 1_33_1 to 1_34_1
File size: 189.0 KB

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

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

Download in other formats: