Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

functional.hpp @ 857

Revision 857, 78.3 KB checked in by igarcia, 18 years ago (diff)

Line
1	//
2	// Copyright (c) 2000-2002
3	// Joerg Walter, Mathias Koch
4	//
5	// Permission to use, copy, modify, distribute and sell this software
6	// and its documentation for any purpose is hereby granted without fee,
7	// provided that the above copyright notice appear in all copies and
8	// that both that copyright notice and this permission notice appear
9	// in supporting documentation. The authors make no representations
10	// about the suitability of this software for any purpose.
11	// It is provided "as is" without express or implied warranty.
12	//
13	// The authors gratefully acknowledge the support of
14	// GeNeSys mbH & Co. KG in producing this work.
15	//
16
17	#ifndef _BOOST_UBLAS_FUNCTIONAL_
18	#define _BOOST_UBLAS_FUNCTIONAL_
19
20	#include <functional>
21
22	#include <boost/numeric/ublas/traits.hpp>
23	#ifdef BOOST_UBLAS_USE_DUFF_DEVICE
24	#include <boost/numeric/ublas/detail/duff.hpp>
25	#endif
26	#ifdef BOOST_UBLAS_USE_SIMD
27	#include <boost/numeric/ublas/detail/raw.hpp>
28	#else
29	namespace boost { namespace numeric { namespace ublas { namespace raw {
30	}}}}
31	#endif
32	#ifdef BOOST_UBLAS_HAVE_BINDINGS
33	#include <boost/numeric/bindings/traits/std_vector.hpp>
34	#include <boost/numeric/bindings/traits/ublas_vector.hpp>
35	#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
36	#include <boost/numeric/bindings/atlas/cblas.hpp>
37	#endif
38
39	#include <boost/numeric/ublas/detail/definitions.hpp>
40
41
42
43	namespace boost { namespace numeric { namespace ublas {
44
45	// Scalar functors
46
47	// Unary
48	template<class T>
49	struct scalar_unary_functor {
50	typedef T value_type;
51	typedef typename type_traits<T>::const_reference argument_type;
52	typedef typename type_traits<T>::value_type result_type;
53	};
54
55	template<class T>
56	struct scalar_identity:
57	public scalar_unary_functor<T> {
58	typedef typename scalar_unary_functor<T>::argument_type argument_type;
59	typedef typename scalar_unary_functor<T>::result_type result_type;
60
61	static BOOST_UBLAS_INLINE
62	result_type apply (argument_type t) {
63	return t;
64	}
65	};
66	template<class T>
67	struct scalar_negate:
68	public scalar_unary_functor<T> {
69	typedef typename scalar_unary_functor<T>::argument_type argument_type;
70	typedef typename scalar_unary_functor<T>::result_type result_type;
71
72	static BOOST_UBLAS_INLINE
73	result_type apply (argument_type t) {
74	return - t;
75	}
76	};
77	template<class T>
78	struct scalar_conj:
79	public scalar_unary_functor<T> {
80	typedef typename scalar_unary_functor<T>::value_type value_type;
81	typedef typename scalar_unary_functor<T>::argument_type argument_type;
82	typedef typename scalar_unary_functor<T>::result_type result_type;
83
84	static BOOST_UBLAS_INLINE
85	result_type apply (argument_type t) {
86	return type_traits<value_type>::conj (t);
87	}
88	};
89
90	// Unary returning real
91	template<class T>
92	struct scalar_real_unary_functor {
93	typedef T value_type;
94	typedef typename type_traits<T>::const_reference argument_type;
95	typedef typename type_traits<T>::real_type result_type;
96	};
97
98	template<class T>
99	struct scalar_real:
100	public scalar_real_unary_functor<T> {
101	typedef typename scalar_real_unary_functor<T>::value_type value_type;
102	typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
103	typedef typename scalar_real_unary_functor<T>::result_type result_type;
104
105	static BOOST_UBLAS_INLINE
106	result_type apply (argument_type t) {
107	return type_traits<value_type>::real (t);
108	}
109	};
110	template<class T>
111	struct scalar_imag:
112	public scalar_real_unary_functor<T> {
113	typedef typename scalar_real_unary_functor<T>::value_type value_type;
114	typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
115	typedef typename scalar_real_unary_functor<T>::result_type result_type;
116
117	static BOOST_UBLAS_INLINE
118	result_type apply (argument_type t) {
119	return type_traits<value_type>::imag (t);
120	}
121	};
122
123	// Binary
124	template<class T1, class T2>
125	struct scalar_binary_functor {
126	typedef typename type_traits<T1>::const_reference argument1_type;
127	typedef typename type_traits<T2>::const_reference argument2_type;
128	typedef typename promote_traits<T1, T2>::promote_type result_type;
129	};
130
131	template<class T1, class T2>
132	struct scalar_plus:
133	public scalar_binary_functor<T1, T2> {
134	typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
135	typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
136	typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
137
138	static BOOST_UBLAS_INLINE
139	result_type apply (argument1_type t1, argument2_type t2) {
140	return t1 + t2;
141	}
142	};
143	template<class T1, class T2>
144	struct scalar_minus:
145	public scalar_binary_functor<T1, T2> {
146	typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
147	typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
148	typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
149
150	static BOOST_UBLAS_INLINE
151	result_type apply (argument1_type t1, argument2_type t2) {
152	return t1 - t2;
153	}
154	};
155	template<class T1, class T2>
156	struct scalar_multiplies:
157	public scalar_binary_functor<T1, T2> {
158	typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
159	typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
160	typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
161
162	static BOOST_UBLAS_INLINE
163	result_type apply (argument1_type t1, argument2_type t2) {
164	return t1 * t2;
165	}
166	};
167	template<class T1, class T2>
168	struct scalar_divides:
169	public scalar_binary_functor<T1, T2> {
170	typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
171	typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
172	typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
173
174	static BOOST_UBLAS_INLINE
175	result_type apply (argument1_type t1, argument2_type t2) {
176	return t1 / t2;
177	}
178	};
179
180	template<class T1, class T2>
181	struct scalar_binary_assign_functor {
182	// ISSUE Remove reference to avoid reference to reference problems
183	typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
184	typedef typename type_traits<T2>::const_reference argument2_type;
185	};
186
187	struct assign_tag {};
188	struct computed_assign_tag {};
189
190	template<class T1, class T2>
191	struct scalar_assign:
192	public scalar_binary_assign_functor<T1, T2> {
193	typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
194	typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
195	#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
196	static const bool computed ;
197	#else
198	static const bool computed = false ;
199	#endif
200
201	static BOOST_UBLAS_INLINE
202	void apply (argument1_type t1, argument2_type t2) {
203	t1 = t2;
204	}
205
206	template<class U1, class U2>
207	struct rebind {
208	typedef scalar_assign<U1, U2> other;
209	};
210	};
211
212	#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
213	template<class T1, class T2>
214	const bool scalar_assign<T1,T2>::computed = false;
215	#endif
216
217	template<class T1, class T2>
218	struct scalar_plus_assign:
219	public scalar_binary_assign_functor<T1, T2> {
220	typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
221	typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
222	#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
223	static const bool computed ;
224	#else
225	static const bool computed = true ;
226	#endif
227
228	static BOOST_UBLAS_INLINE
229	void apply (argument1_type t1, argument2_type t2) {
230	t1 += t2;
231	}
232
233	template<class U1, class U2>
234	struct rebind {
235	typedef scalar_plus_assign<U1, U2> other;
236	};
237	};
238
239	#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
240	template<class T1, class T2>
241	const bool scalar_plus_assign<T1,T2>::computed = true;
242	#endif
243
244	template<class T1, class T2>
245	struct scalar_minus_assign:
246	public scalar_binary_assign_functor<T1, T2> {
247	typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
248	typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
249	#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
250	static const bool computed ;
251	#else
252	static const bool computed = true ;
253	#endif
254
255	static BOOST_UBLAS_INLINE
256	void apply (argument1_type t1, argument2_type t2) {
257	t1 -= t2;
258	}
259
260	template<class U1, class U2>
261	struct rebind {
262	typedef scalar_minus_assign<U1, U2> other;
263	};
264	};
265
266	#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
267	template<class T1, class T2>
268	const bool scalar_minus_assign<T1,T2>::computed = true;
269	#endif
270
271	template<class T1, class T2>
272	struct scalar_multiplies_assign:
273	public scalar_binary_assign_functor<T1, T2> {
274	typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
275	typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
276	static const bool computed = true;
277
278	static BOOST_UBLAS_INLINE
279	void apply (argument1_type t1, argument2_type t2) {
280	t1 *= t2;
281	}
282
283	template<class U1, class U2>
284	struct rebind {
285	typedef scalar_multiplies_assign<U1, U2> other;
286	};
287	};
288	template<class T1, class T2>
289	struct scalar_divides_assign:
290	public scalar_binary_assign_functor<T1, T2> {
291	typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
292	typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
293	static const bool computed ;
294
295	static BOOST_UBLAS_INLINE
296	void apply (argument1_type t1, argument2_type t2) {
297	t1 /= t2;
298	}
299
300	template<class U1, class U2>
301	struct rebind {
302	typedef scalar_divides_assign<U1, U2> other;
303	};
304	};
305	template<class T1, class T2>
306	const bool scalar_divides_assign<T1,T2>::computed = true;
307
308	template<class T1, class T2>
309	struct scalar_binary_swap_functor {
310	typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
311	typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type;
312	};
313
314	template<class T1, class T2>
315	struct scalar_swap:
316	public scalar_binary_swap_functor<T1, T2> {
317	typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type;
318	typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type;
319
320	static BOOST_UBLAS_INLINE
321	void apply (argument1_type t1, argument2_type t2) {
322	std::swap (t1, t2);
323	}
324
325	template<class U1, class U2>
326	struct rebind {
327	typedef scalar_swap<U1, U2> other;
328	};
329	};
330
331	// Vector functors
332
333	// Unary returning scalar
334	template<class T>
335	struct vector_scalar_unary_functor {
336	typedef std::size_t size_type;
337	typedef std::ptrdiff_t difference_type;
338	typedef T value_type;
339	typedef T result_type;
340	};
341
342	template<class T>
343	struct vector_sum:
344	public vector_scalar_unary_functor<T> {
345	typedef typename vector_scalar_unary_functor<T>::size_type size_type;
346	typedef typename vector_scalar_unary_functor<T>::difference_type difference_type;
347	typedef typename vector_scalar_unary_functor<T>::value_type value_type;
348	typedef typename vector_scalar_unary_functor<T>::result_type result_type;
349
350	template<class E>
351	static BOOST_UBLAS_INLINE
352	result_type apply (const vector_expression<E> &e) {
353	result_type t = result_type (0);
354	size_type size (e ().size ());
355	for (size_type i = 0; i < size; ++ i)
356	t += e () (i);
357	return t;
358	}
359	// Dense case
360	template<class I>
361	static BOOST_UBLAS_INLINE
362	result_type apply (difference_type size, I it) {
363	result_type t = result_type (0);
364	while (-- size >= 0)
365	t += *it, ++ it;
366	return t;
367	}
368	// Sparse case
369	template<class I>
370	static BOOST_UBLAS_INLINE
371	result_type apply (I it, const I &it_end) {
372	result_type t = result_type (0);
373	while (it != it_end)
374	t += *it, ++ it;
375	return t;
376	}
377	};
378
379	// Unary returning real scalar
380	template<class T>
381	struct vector_scalar_real_unary_functor {
382	typedef std::size_t size_type;
383	typedef std::ptrdiff_t difference_type;
384	typedef T value_type;
385	typedef typename type_traits<T>::real_type real_type;
386	typedef real_type result_type;
387	};
388
389	template<class T>
390	struct vector_norm_1:
391	public vector_scalar_real_unary_functor<T> {
392	typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
393	typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
394	typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
395	typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
396	typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
397
398	template<class E>
399	static BOOST_UBLAS_INLINE
400	result_type apply (const vector_expression<E> &e) {
401	real_type t = real_type ();
402	size_type size (e ().size ());
403	for (size_type i = 0; i < size; ++ i) {
404	real_type u (type_traits<value_type>::norm_1 (e () (i)));
405	t += u;
406	}
407	return t;
408	}
409	// Dense case
410	template<class I>
411	static BOOST_UBLAS_INLINE
412	result_type apply (difference_type size, I it) {
413	real_type t = real_type ();
414	while (-- size >= 0) {
415	real_type u (type_traits<value_type>::norm_1 (*it));
416	t += u;
417	++ it;
418	}
419	return t;
420	}
421	// Sparse case
422	template<class I>
423	static BOOST_UBLAS_INLINE
424	result_type apply (I it, const I &it_end) {
425	real_type t = real_type ();
426	while (it != it_end) {
427	real_type u (type_traits<value_type>::norm_1 (*it));
428	t += u;
429	++ it;
430	}
431	return t;
432	}
433	};
434	template<class T>
435	struct vector_norm_2:
436	public vector_scalar_real_unary_functor<T> {
437	typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
438	typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
439	typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
440	typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
441	typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
442
443	template<class E>
444	static BOOST_UBLAS_INLINE
445	result_type apply (const vector_expression<E> &e) {
446	#ifndef BOOST_UBLAS_SCALED_NORM
447	real_type t = real_type ();
448	size_type size (e ().size ());
449	for (size_type i = 0; i < size; ++ i) {
450	real_type u (type_traits<value_type>::norm_2 (e () (i)));
451	t += u * u;
452	}
453	return type_traits<real_type>::sqrt (t);
454	#else
455	real_type scale = real_type ();
456	real_type sum_squares (1);
457	size_type size (e ().size ());
458	for (size_type i = 0; i < size; ++ i) {
459	real_type u (type_traits<value_type>::norm_2 (e () (i)));
460	if (scale < u) {
461	real_type v (scale / u);
462	sum_squares = sum_squares * v * v + real_type (1);
463	scale = u;
464	} else {
465	real_type v (u / scale);
466	sum_squares += v * v;
467	}
468	}
469	return scale * type_traits<real_type>::sqrt (sum_squares);
470	#endif
471	}
472	// Dense case
473	template<class I>
474	static BOOST_UBLAS_INLINE
475	result_type apply (difference_type size, I it) {
476	#ifndef BOOST_UBLAS_SCALED_NORM
477	real_type t = real_type ();
478	while (-- size >= 0) {
479	real_type u (type_traits<value_type>::norm_2 (*it));
480	t += u * u;
481	++ it;
482	}
483	return type_traits<real_type>::sqrt (t);
484	#else
485	real_type scale = real_type ();
486	real_type sum_squares (1);
487	while (-- size >= 0) {
488	real_type u (type_traits<value_type>::norm_2 (*it));
489	if (scale < u) {
490	real_type v (scale / u);
491	sum_squares = sum_squares * v * v + real_type (1);
492	scale = u;
493	} else {
494	real_type v (u / scale);
495	sum_squares += v * v;
496	}
497	++ it;
498	}
499	return scale * type_traits<real_type>::sqrt (sum_squares);
500	#endif
501	}
502	// Sparse case
503	template<class I>
504	static BOOST_UBLAS_INLINE
505	result_type apply (I it, const I &it_end) {
506	#ifndef BOOST_UBLAS_SCALED_NORM
507	real_type t = real_type ();
508	while (it != it_end) {
509	real_type u (type_traits<value_type>::norm_2 (*it));
510	t += u * u;
511	++ it;
512	}
513	return type_traits<real_type>::sqrt (t);
514	#else
515	real_type scale = real_type ();
516	real_type sum_squares (1);
517	while (it != it_end) {
518	real_type u (type_traits<value_type>::norm_2 (*it));
519	if (scale < u) {
520	real_type v (scale / u);
521	sum_squares = sum_squares * v * v + real_type (1);
522	scale = u;
523	} else {
524	real_type v (u / scale);
525	sum_squares += v * v;
526	}
527	++ it;
528	}
529	return scale * type_traits<real_type>::sqrt (sum_squares);
530	#endif
531	}
532	};
533	template<class T>
534	struct vector_norm_inf:
535	public vector_scalar_real_unary_functor<T> {
536	typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
537	typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
538	typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
539	typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
540	typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
541
542	template<class E>
543	static BOOST_UBLAS_INLINE
544	result_type apply (const vector_expression<E> &e) {
545	real_type t = real_type ();
546	size_type size (e ().size ());
547	for (size_type i = 0; i < size; ++ i) {
548	real_type u (type_traits<value_type>::norm_inf (e () (i)));
549	if (u > t)
550	t = u;
551	}
552	return t;
553	}
554	// Dense case
555	template<class I>
556	static BOOST_UBLAS_INLINE
557	result_type apply (difference_type size, I it) {
558	real_type t = real_type ();
559	while (-- size >= 0) {
560	real_type u (type_traits<value_type>::norm_inf (*it));
561	if (u > t)
562	t = u;
563	++ it;
564	}
565	return t;
566	}
567	// Sparse case
568	template<class I>
569	static BOOST_UBLAS_INLINE
570	result_type apply (I it, const I &it_end) {
571	real_type t = real_type ();
572	while (it != it_end) {
573	real_type u (type_traits<value_type>::norm_inf (*it));
574	if (u > t)
575	t = u;
576	++ it;
577	}
578	return t;
579	}
580	};
581
582	// Unary returning index
583	template<class T>
584	struct vector_scalar_index_unary_functor {
585	typedef std::size_t size_type;
586	typedef std::ptrdiff_t difference_type;
587	typedef T value_type;
588	typedef typename type_traits<T>::real_type real_type;
589	typedef size_type result_type;
590	};
591
592	template<class T>
593	struct vector_index_norm_inf:
594	public vector_scalar_index_unary_functor<T> {
595	typedef typename vector_scalar_index_unary_functor<T>::size_type size_type;
596	typedef typename vector_scalar_index_unary_functor<T>::difference_type difference_type;
597	typedef typename vector_scalar_index_unary_functor<T>::value_type value_type;
598	typedef typename vector_scalar_index_unary_functor<T>::real_type real_type;
599	typedef typename vector_scalar_index_unary_functor<T>::result_type result_type;
600
601	template<class E>
602	static BOOST_UBLAS_INLINE
603	result_type apply (const vector_expression<E> &e) {
604	// ISSUE For CBLAS compatibility return 0 index in empty case
605	result_type i_norm_inf (0);
606	real_type t = real_type ();
607	size_type size (e ().size ());
608	for (size_type i = 0; i < size; ++ i) {
609	real_type u (type_traits<value_type>::norm_inf (e () (i)));
610	if (u > t) {
611	i_norm_inf = i;
612	t = u;
613	}
614	}
615	return i_norm_inf;
616	}
617	// Dense case
618	template<class I>
619	static BOOST_UBLAS_INLINE
620	result_type apply (difference_type size, I it) {
621	// ISSUE For CBLAS compatibility return 0 index in empty case
622	result_type i_norm_inf (0);
623	real_type t = real_type ();
624	while (-- size >= 0) {
625	real_type u (type_traits<value_type>::norm_inf (*it));
626	if (u > t) {
627	i_norm_inf = it.index ();
628	t = u;
629	}
630	++ it;
631	}
632	return i_norm_inf;
633	}
634	// Sparse case
635	template<class I>
636	static BOOST_UBLAS_INLINE
637	result_type apply (I it, const I &it_end) {
638	// ISSUE For CBLAS compatibility return 0 index in empty case
639	result_type i_norm_inf (0);
640	real_type t = real_type ();
641	while (it != it_end) {
642	real_type u (type_traits<value_type>::norm_inf (*it));
643	if (u > t) {
644	i_norm_inf = it.index ();
645	t = u;
646	}
647	++ it;
648	}
649	return i_norm_inf;
650	}
651	};
652
653	// Binary returning scalar
654	template<class T1, class T2, class TR>
655	struct vector_scalar_binary_functor {
656	typedef std::size_t size_type;
657	typedef std::ptrdiff_t difference_type;
658	typedef TR value_type;
659	typedef TR result_type;
660	};
661
662	template<class T1, class T2, class TR>
663	struct vector_inner_prod:
664	public vector_scalar_binary_functor<T1, T2, TR> {
665	typedef typename vector_scalar_binary_functor<T1, T2, TR>::size_type size_type ;
666	typedef typename vector_scalar_binary_functor<T1, T2, TR>::difference_type difference_type;
667	typedef typename vector_scalar_binary_functor<T1, T2, TR>::value_type value_type;
668	typedef typename vector_scalar_binary_functor<T1, T2, TR>::result_type result_type;
669
670	template<class C1, class C2>
671	static BOOST_UBLAS_INLINE
672	result_type apply (const vector_container<C1> &c1,
673	const vector_container<C2> &c2) {
674	#ifdef BOOST_UBLAS_USE_SIMD
675	using namespace raw;
676	size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
677	const T1 *data1 = data_const (c1 ());
678	const T2 *data2 = data_const (c2 ());
679	size_type s1 = stride (c1 ());
680	size_type s2 = stride (c2 ());
681	result_type t = result_type (0);
682	if (s1 == 1 && s2 == 1) {
683	for (size_type i = 0; i < size; ++ i)
684	t += data1 [i] * data2 [i];
685	} else if (s2 == 1) {
686	for (size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
687	t += data1 [i1] * data2 [i];
688	} else if (s1 == 1) {
689	for (size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
690	t += data1 [i] * data2 [i2];
691	} else {
692	for (size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
693	t += data1 [i1] * data2 [i2];
694	}
695	return t;
696	#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
697	return boost::numeric::bindings::atlas::dot (c1 (), c2 ());
698	#else
699	return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2));
700	#endif
701	}
702	template<class E1, class E2>
703	static BOOST_UBLAS_INLINE
704	result_type apply (const vector_expression<E1> &e1,
705	const vector_expression<E2> &e2) {
706	size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
707	result_type t = result_type (0);
708	#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
709	for (size_type i = 0; i < size; ++ i)
710	t += e1 () (i) * e2 () (i);
711	#else
712	size_type i (0);
713	DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));
714	#endif
715	return t;
716	}
717	// Dense case
718	template<class I1, class I2>
719	static BOOST_UBLAS_INLINE
720	result_type apply (difference_type size, I1 it1, I2 it2) {
721	result_type t = result_type (0);
722	#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
723	while (-- size >= 0)
724	t += it1 *it2, ++ it1, ++ it2;
725	#else
726	DD (size, 4, r, (t += it1 *it2, ++ it1, ++ it2));
727	#endif
728	return t;
729	}
730	// Packed case
731	template<class I1, class I2>
732	static BOOST_UBLAS_INLINE
733	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
734	result_type t = result_type (0);
735	difference_type it1_size (it1_end - it1);
736	difference_type it2_size (it2_end - it2);
737	difference_type diff (0);
738	if (it1_size > 0 && it2_size > 0)
739	diff = it2.index () - it1.index ();
740	if (diff != 0) {
741	difference_type size = (std::min) (diff, it1_size);
742	if (size > 0) {
743	it1 += size;
744	it1_size -= size;
745	diff -= size;
746	}
747	size = (std::min) (- diff, it2_size);
748	if (size > 0) {
749	it2 += size;
750	it2_size -= size;
751	diff += size;
752	}
753	}
754	difference_type size ((std::min) (it1_size, it2_size));
755	while (-- size >= 0)
756	t += it1 *it2, ++ it1, ++ it2;
757	return t;
758	}
759	// Sparse case
760	template<class I1, class I2>
761	static BOOST_UBLAS_INLINE
762	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
763	result_type t = result_type (0);
764	if (it1 != it1_end && it2 != it2_end) {
765	size_type it1_index = it1.index (), it2_index = it2.index ();
766	while (true) {
767	difference_type compare = it1_index - it2_index;
768	if (compare == 0) {
769	t += it1 *it2, ++ it1, ++ it2;
770	if (it1 != it1_end && it2 != it2_end) {
771	it1_index = it1.index ();
772	it2_index = it2.index ();
773	} else
774	break;
775	} else if (compare < 0) {
776	increment (it1, it1_end, - compare);
777	if (it1 != it1_end)
778	it1_index = it1.index ();
779	else
780	break;
781	} else if (compare > 0) {
782	increment (it2, it2_end, compare);
783	if (it2 != it2_end)
784	it2_index = it2.index ();
785	else
786	break;
787	}
788	}
789	}
790	return t;
791	}
792	};
793
794	// Matrix functors
795
796	// Binary returning vector
797	template<class T1, class T2, class TR>
798	struct matrix_vector_binary_functor {
799	typedef std::size_t size_type;
800	typedef std::ptrdiff_t difference_type;
801	typedef TR value_type;
802	typedef TR result_type;
803	};
804
805	template<class T1, class T2, class TR>
806	struct matrix_vector_prod1:
807	public matrix_vector_binary_functor<T1, T2, TR> {
808	typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;
809	typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;
810	typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;
811	typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;
812
813	template<class C1, class C2>
814	static BOOST_UBLAS_INLINE
815	result_type apply (const matrix_container<C1> &c1,
816	const vector_container<C2> &c2,
817	size_type i) {
818	#ifdef BOOST_UBLAS_USE_SIMD
819	using namespace raw;
820	size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ());
821	const T1 data1 = data_const (c1 ()) + i stride1 (c1 ());
822	const T2 *data2 = data_const (c2 ());
823	size_type s1 = stride2 (c1 ());
824	size_type s2 = stride (c2 ());
825	result_type t = result_type (0);
826	if (s1 == 1 && s2 == 1) {
827	for (size_type j = 0; j < size; ++ j)
828	t += data1 [j] * data2 [j];
829	} else if (s2 == 1) {
830	for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
831	t += data1 [j1] * data2 [j];
832	} else if (s1 == 1) {
833	for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
834	t += data1 [j] * data2 [j2];
835	} else {
836	for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
837	t += data1 [j1] * data2 [j2];
838	}
839	return t;
840	#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
841	return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ());
842	#else
843	return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i));
844	#endif
845	}
846	template<class E1, class E2>
847	static BOOST_UBLAS_INLINE
848	result_type apply (const matrix_expression<E1> &e1,
849	const vector_expression<E2> &e2,
850	size_type i) {
851	size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
852	result_type t = result_type (0);
853	#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
854	for (size_type j = 0; j < size; ++ j)
855	t += e1 () (i, j) * e2 () (j);
856	#else
857	size_type j (0);
858	DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j));
859	#endif
860	return t;
861	}
862	// Dense case
863	template<class I1, class I2>
864	static BOOST_UBLAS_INLINE
865	result_type apply (difference_type size, I1 it1, I2 it2) {
866	result_type t = result_type (0);
867	#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
868	while (-- size >= 0)
869	t += it1 *it2, ++ it1, ++ it2;
870	#else
871	DD (size, 4, r, (t += it1 *it2, ++ it1, ++ it2));
872	#endif
873	return t;
874	}
875	// Packed case
876	template<class I1, class I2>
877	static BOOST_UBLAS_INLINE
878	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
879	result_type t = result_type (0);
880	difference_type it1_size (it1_end - it1);
881	difference_type it2_size (it2_end - it2);
882	difference_type diff (0);
883	if (it1_size > 0 && it2_size > 0)
884	diff = it2.index () - it1.index2 ();
885	if (diff != 0) {
886	difference_type size = (std::min) (diff, it1_size);
887	if (size > 0) {
888	it1 += size;
889	it1_size -= size;
890	diff -= size;
891	}
892	size = (std::min) (- diff, it2_size);
893	if (size > 0) {
894	it2 += size;
895	it2_size -= size;
896	diff += size;
897	}
898	}
899	difference_type size ((std::min) (it1_size, it2_size));
900	while (-- size >= 0)
901	t += it1 *it2, ++ it1, ++ it2;
902	return t;
903	}
904	// Sparse case
905	template<class I1, class I2>
906	static BOOST_UBLAS_INLINE
907	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
908	sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
909	result_type t = result_type (0);
910	if (it1 != it1_end && it2 != it2_end) {
911	size_type it1_index = it1.index2 (), it2_index = it2.index ();
912	while (true) {
913	difference_type compare = it1_index - it2_index;
914	if (compare == 0) {
915	t += it1 *it2, ++ it1, ++ it2;
916	if (it1 != it1_end && it2 != it2_end) {
917	it1_index = it1.index2 ();
918	it2_index = it2.index ();
919	} else
920	break;
921	} else if (compare < 0) {
922	increment (it1, it1_end, - compare);
923	if (it1 != it1_end)
924	it1_index = it1.index2 ();
925	else
926	break;
927	} else if (compare > 0) {
928	increment (it2, it2_end, compare);
929	if (it2 != it2_end)
930	it2_index = it2.index ();
931	else
932	break;
933	}
934	}
935	}
936	return t;
937	}
938	// Sparse packed case
939	template<class I1, class I2>
940	static BOOST_UBLAS_INLINE
941	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
942	sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
943	result_type t = result_type (0);
944	while (it1 != it1_end) {
945	t += it1 it2 () (it1.index2 ());
946	++ it1;
947	}
948	return t;
949	}
950	// Packed sparse case
951	template<class I1, class I2>
952	static BOOST_UBLAS_INLINE
953	result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
954	packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
955	result_type t = result_type (0);
956	while (it2 != it2_end) {
957	t += it1 () (it1.index1 (), it2.index ()) * *it2;
958	++ it2;
959	}
960	return t;
961	}
962	// Another dispatcher
963	template<class I1, class I2>
964	static BOOST_UBLAS_INLINE
965	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
966	sparse_bidirectional_iterator_tag) {
967	typedef typename I1::iterator_category iterator1_category;
968	typedef typename I2::iterator_category iterator2_category;
969	return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
970	}
971	};
972
973	template<class T1, class T2, class TR>
974	struct matrix_vector_prod2:
975	public matrix_vector_binary_functor<T1, T2, TR> {
976	typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;
977	typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;
978	typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;
979	typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;
980
981	template<class C1, class C2>
982	static BOOST_UBLAS_INLINE
983	result_type apply (const vector_container<C1> &c1,
984	const matrix_container<C2> &c2,
985	size_type i) {
986	#ifdef BOOST_UBLAS_USE_SIMD
987	using namespace raw;
988	size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ());
989	const T1 *data1 = data_const (c1 ());
990	const T2 data2 = data_const (c2 ()) + i stride2 (c2 ());
991	size_type s1 = stride (c1 ());
992	size_type s2 = stride1 (c2 ());
993	result_type t = result_type (0);
994	if (s1 == 1 && s2 == 1) {
995	for (size_type j = 0; j < size; ++ j)
996	t += data1 [j] * data2 [j];
997	} else if (s2 == 1) {
998	for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
999	t += data1 [j1] * data2 [j];
1000	} else if (s1 == 1) {
1001	for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
1002	t += data1 [j] * data2 [j2];
1003	} else {
1004	for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
1005	t += data1 [j1] * data2 [j2];
1006	}
1007	return t;
1008	#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1009	return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i));
1010	#else
1011	return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
1012	#endif
1013	}
1014	template<class E1, class E2>
1015	static BOOST_UBLAS_INLINE
1016	result_type apply (const vector_expression<E1> &e1,
1017	const matrix_expression<E2> &e2,
1018	size_type i) {
1019	size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
1020	result_type t = result_type (0);
1021	#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1022	for (size_type j = 0; j < size; ++ j)
1023	t += e1 () (j) * e2 () (j, i);
1024	#else
1025	size_type j (0);
1026	DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j));
1027	#endif
1028	return t;
1029	}
1030	// Dense case
1031	template<class I1, class I2>
1032	static BOOST_UBLAS_INLINE
1033	result_type apply (difference_type size, I1 it1, I2 it2) {
1034	result_type t = result_type (0);
1035	#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1036	while (-- size >= 0)
1037	t += it1 *it2, ++ it1, ++ it2;
1038	#else
1039	DD (size, 4, r, (t += it1 *it2, ++ it1, ++ it2));
1040	#endif
1041	return t;
1042	}
1043	// Packed case
1044	template<class I1, class I2>
1045	static BOOST_UBLAS_INLINE
1046	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
1047	result_type t = result_type (0);
1048	difference_type it1_size (it1_end - it1);
1049	difference_type it2_size (it2_end - it2);
1050	difference_type diff (0);
1051	if (it1_size > 0 && it2_size > 0)
1052	diff = it2.index1 () - it1.index ();
1053	if (diff != 0) {
1054	difference_type size = (std::min) (diff, it1_size);
1055	if (size > 0) {
1056	it1 += size;
1057	it1_size -= size;
1058	diff -= size;
1059	}
1060	size = (std::min) (- diff, it2_size);
1061	if (size > 0) {
1062	it2 += size;
1063	it2_size -= size;
1064	diff += size;
1065	}
1066	}
1067	difference_type size ((std::min) (it1_size, it2_size));
1068	while (-- size >= 0)
1069	t += it1 *it2, ++ it1, ++ it2;
1070	return t;
1071	}
1072	// Sparse case
1073	template<class I1, class I2>
1074	static BOOST_UBLAS_INLINE
1075	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1076	sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
1077	result_type t = result_type (0);
1078	if (it1 != it1_end && it2 != it2_end) {
1079	size_type it1_index = it1.index (), it2_index = it2.index1 ();
1080	while (true) {
1081	difference_type compare = it1_index - it2_index;
1082	if (compare == 0) {
1083	t += it1 *it2, ++ it1, ++ it2;
1084	if (it1 != it1_end && it2 != it2_end) {
1085	it1_index = it1.index ();
1086	it2_index = it2.index1 ();
1087	} else
1088	break;
1089	} else if (compare < 0) {
1090	increment (it1, it1_end, - compare);
1091	if (it1 != it1_end)
1092	it1_index = it1.index ();
1093	else
1094	break;
1095	} else if (compare > 0) {
1096	increment (it2, it2_end, compare);
1097	if (it2 != it2_end)
1098	it2_index = it2.index1 ();
1099	else
1100	break;
1101	}
1102	}
1103	}
1104	return t;
1105	}
1106	// Packed sparse case
1107	template<class I1, class I2>
1108	static BOOST_UBLAS_INLINE
1109	result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
1110	packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
1111	result_type t = result_type (0);
1112	while (it2 != it2_end) {
1113	t += it1 () (it2.index1 ()) * *it2;
1114	++ it2;
1115	}
1116	return t;
1117	}
1118	// Sparse packed case
1119	template<class I1, class I2>
1120	static BOOST_UBLAS_INLINE
1121	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
1122	sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
1123	result_type t = result_type (0);
1124	while (it1 != it1_end) {
1125	t += it1 it2 () (it1.index (), it2.index2 ());
1126	++ it1;
1127	}
1128	return t;
1129	}
1130	// Another dispatcher
1131	template<class I1, class I2>
1132	static BOOST_UBLAS_INLINE
1133	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1134	sparse_bidirectional_iterator_tag) {
1135	typedef typename I1::iterator_category iterator1_category;
1136	typedef typename I2::iterator_category iterator2_category;
1137	return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
1138	}
1139	};
1140
1141	// Binary returning matrix
1142	template<class T1, class T2, class TR>
1143	struct matrix_matrix_binary_functor {
1144	typedef std::size_t size_type;
1145	typedef std::ptrdiff_t difference_type;
1146	typedef TR value_type;
1147	typedef TR result_type;
1148	};
1149
1150	template<class T1, class T2, class TR>
1151	struct matrix_matrix_prod:
1152	public matrix_matrix_binary_functor<T1, T2, TR> {
1153	typedef typename matrix_matrix_binary_functor<T1, T2, TR>::size_type size_type;
1154	typedef typename matrix_matrix_binary_functor<T1, T2, TR>::difference_type difference_type;
1155	typedef typename matrix_matrix_binary_functor<T1, T2, TR>::value_type value_type;
1156	typedef typename matrix_matrix_binary_functor<T1, T2, TR>::result_type result_type;
1157
1158	template<class C1, class C2>
1159	static BOOST_UBLAS_INLINE
1160	result_type apply (const matrix_container<C1> &c1,
1161	const matrix_container<C2> &c2,
1162	size_type i, size_type j) {
1163	#ifdef BOOST_UBLAS_USE_SIMD
1164	using namespace raw;
1165	size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ());
1166	const T1 data1 = data_const (c1 ()) + i stride1 (c1 ());
1167	const T2 data2 = data_const (c2 ()) + j stride2 (c2 ());
1168	size_type s1 = stride2 (c1 ());
1169	size_type s2 = stride1 (c2 ());
1170	result_type t = result_type (0);
1171	if (s1 == 1 && s2 == 1) {
1172	for (size_type k = 0; k < size; ++ k)
1173	t += data1 [k] * data2 [k];
1174	} else if (s2 == 1) {
1175	for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1)
1176	t += data1 [k1] * data2 [k];
1177	} else if (s1 == 1) {
1178	for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2)
1179	t += data1 [k] * data2 [k2];
1180	} else {
1181	for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2)
1182	t += data1 [k1] * data2 [k2];
1183	}
1184	return t;
1185	#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1186	return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j));
1187	#else
1188	return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
1189	#endif
1190	}
1191	template<class E1, class E2>
1192	static BOOST_UBLAS_INLINE
1193	result_type apply (const matrix_expression<E1> &e1,
1194	const matrix_expression<E2> &e2,
1195	size_type i, size_type j) {
1196	size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
1197	result_type t = result_type (0);
1198	#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1199	for (size_type k = 0; k < size; ++ k)
1200	t += e1 () (i, k) * e2 () (k, j);
1201	#else
1202	size_type k (0);
1203	DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k));
1204	#endif
1205	return t;
1206	}
1207	// Dense case
1208	template<class I1, class I2>
1209	static BOOST_UBLAS_INLINE
1210	result_type apply (difference_type size, I1 it1, I2 it2) {
1211	result_type t = result_type (0);
1212	#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1213	while (-- size >= 0)
1214	t += it1 *it2, ++ it1, ++ it2;
1215	#else
1216	DD (size, 4, r, (t += it1 *it2, ++ it1, ++ it2));
1217	#endif
1218	return t;
1219	}
1220	// Packed case
1221	template<class I1, class I2>
1222	static BOOST_UBLAS_INLINE
1223	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) {
1224	result_type t = result_type (0);
1225	difference_type it1_size (it1_end - it1);
1226	difference_type it2_size (it2_end - it2);
1227	difference_type diff (0);
1228	if (it1_size > 0 && it2_size > 0)
1229	diff = it2.index1 () - it1.index2 ();
1230	if (diff != 0) {
1231	difference_type size = (std::min) (diff, it1_size);
1232	if (size > 0) {
1233	it1 += size;
1234	it1_size -= size;
1235	diff -= size;
1236	}
1237	size = (std::min) (- diff, it2_size);
1238	if (size > 0) {
1239	it2 += size;
1240	it2_size -= size;
1241	diff += size;
1242	}
1243	}
1244	difference_type size ((std::min) (it1_size, it2_size));
1245	while (-- size >= 0)
1246	t += it1 *it2, ++ it1, ++ it2;
1247	return t;
1248	}
1249	// Sparse case
1250	template<class I1, class I2>
1251	static BOOST_UBLAS_INLINE
1252	result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
1253	result_type t = result_type (0);
1254	if (it1 != it1_end && it2 != it2_end) {
1255	size_type it1_index = it1.index2 (), it2_index = it2.index1 ();
1256	while (true) {
1257	difference_type compare = it1_index - it2_index;
1258	if (compare == 0) {
1259	t += it1 *it2, ++ it1, ++ it2;
1260	if (it1 != it1_end && it2 != it2_end) {
1261	it1_index = it1.index2 ();
1262	it2_index = it2.index1 ();
1263	} else
1264	break;
1265	} else if (compare < 0) {
1266	increment (it1, it1_end, - compare);
1267	if (it1 != it1_end)
1268	it1_index = it1.index2 ();
1269	else
1270	break;
1271	} else if (compare > 0) {
1272	increment (it2, it2_end, compare);
1273	if (it2 != it2_end)
1274	it2_index = it2.index1 ();
1275	else
1276	break;
1277	}
1278	}
1279	}
1280	return t;
1281	}
1282	};
1283
1284	// Unary returning scalar norm
1285	template<class T>
1286	struct matrix_scalar_real_unary_functor {
1287	typedef std::size_t size_type;
1288	typedef std::ptrdiff_t difference_type;
1289	typedef T value_type;
1290	typedef typename type_traits<T>::real_type real_type;
1291	typedef real_type result_type;
1292	};
1293
1294	template<class T>
1295	struct matrix_norm_1:
1296	public matrix_scalar_real_unary_functor<T> {
1297	typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
1298	typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
1299	typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
1300	typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
1301	typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
1302
1303	template<class E>
1304	static BOOST_UBLAS_INLINE
1305	result_type apply (const matrix_expression<E> &e) {
1306	real_type t = real_type ();
1307	size_type size2 (e ().size2 ());
1308	for (size_type j = 0; j < size2; ++ j) {
1309	real_type u = real_type ();
1310	size_type size1 (e ().size1 ());
1311	for (size_type i = 0; i < size1; ++ i) {
1312	real_type v (type_traits<value_type>::norm_1 (e () (i, j)));
1313	u += v;
1314	}
1315	if (u > t)
1316	t = u;
1317	}
1318	return t;
1319	}
1320	};
1321	template<class T>
1322	struct matrix_norm_frobenius:
1323	public matrix_scalar_real_unary_functor<T> {
1324	typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
1325	typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
1326	typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
1327	typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
1328	typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
1329
1330	template<class E>
1331	static BOOST_UBLAS_INLINE
1332	result_type apply (const matrix_expression<E> &e) {
1333	real_type t = real_type ();
1334	size_type size1 (e ().size1 ());
1335	for (size_type i = 0; i < size1; ++ i) {
1336	size_type size2 (e ().size2 ());
1337	for (size_type j = 0; j < size2; ++ j) {
1338	real_type u (type_traits<value_type>::norm_2 (e () (i, j)));
1339	t += u * u;
1340	}
1341	}
1342	return type_traits<real_type>::sqrt (t);
1343	}
1344	};
1345	template<class T>
1346	struct matrix_norm_inf:
1347	public matrix_scalar_real_unary_functor<T> {
1348	typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
1349	typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
1350	typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
1351	typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
1352	typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
1353
1354	template<class E>
1355	static BOOST_UBLAS_INLINE
1356	result_type apply (const matrix_expression<E> &e) {
1357	real_type t = real_type ();
1358	size_type size1 (e ().size1 ());
1359	for (size_type i = 0; i < size1; ++ i) {
1360	real_type u = real_type ();
1361	size_type size2 (e ().size2 ());
1362	for (size_type j = 0; j < size2; ++ j) {
1363	real_type v (type_traits<value_type>::norm_inf (e () (i, j)));
1364	u += v;
1365	}
1366	if (u > t)
1367	t = u;
1368	}
1369	return t;
1370	}
1371	};
1372
1373	// This functor computes the address translation
1374	// matrix [i] [j] -> storage [i * size2 + j]
1375	template <class Z, class D>
1376	struct basic_row_major {
1377	typedef Z size_type;
1378	typedef D difference_type;
1379	typedef row_major_tag orientation_category;
1380
1381	static
1382	BOOST_UBLAS_INLINE
1383	size_type storage_size (size_type size1, size_type size2) {
1384	// Guard against size_type overflow
1385	BOOST_UBLAS_CHECK (size2 == 0 \|\| size1 <= (std::numeric_limits<size_type>::max) () / size2, bad_size ());
1386	return size1 * size2;
1387	}
1388
1389	// Indexing
1390	static
1391	BOOST_UBLAS_INLINE
1392	size_type element (size_type i, size_type size1, size_type j, size_type size2) {
1393	BOOST_UBLAS_CHECK (i < size1, bad_index ());
1394	BOOST_UBLAS_CHECK (j < size2, bad_index ());
1395	detail::ignore_unused_variable_warning(size1);
1396	// Guard against size_type overflow
1397	BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size2, bad_index ());
1398	return i * size2 + j;
1399	}
1400	static
1401	BOOST_UBLAS_INLINE
1402	size_type address (size_type i, size_type size1, size_type j, size_type size2) {
1403	BOOST_UBLAS_CHECK (i <= size1, bad_index ());
1404	BOOST_UBLAS_CHECK (j <= size2, bad_index ());
1405	// Guard against size_type overflow - address may be size2 past end of storage
1406	BOOST_UBLAS_CHECK (size2 == 0 \|\| i <= ((std::numeric_limits<size_type>::max) () - j) / size2, bad_index ());
1407	detail::ignore_unused_variable_warning(size1);
1408	return i * size2 + j;
1409	}
1410	static
1411	BOOST_UBLAS_INLINE
1412	difference_type distance1 (difference_type k, size_type /* size1 */, size_type size2) {
1413	return size2 != 0 ? k / size2 : 0;
1414	}
1415	static
1416	BOOST_UBLAS_INLINE
1417	difference_type distance2 (difference_type k, size_type /* size1 /, size_type / size2 */) {
1418	return k;
1419	}
1420	static
1421	BOOST_UBLAS_INLINE
1422	size_type index1 (difference_type k, size_type /* size1 */, size_type size2) {
1423	return size2 != 0 ? k / size2 : 0;
1424	}
1425	static
1426	BOOST_UBLAS_INLINE
1427	size_type index2 (difference_type k, size_type /* size1 */, size_type size2) {
1428	return size2 != 0 ? k % size2 : 0;
1429	}
1430	static
1431	BOOST_UBLAS_INLINE
1432	bool fast1 () {
1433	return false;
1434	}
1435	static
1436	BOOST_UBLAS_INLINE
1437	size_type one1 (size_type /* size1 */, size_type size2) {
1438	return size2;
1439	}
1440	static
1441	BOOST_UBLAS_INLINE
1442	bool fast2 () {
1443	return true;
1444	}
1445	static
1446	BOOST_UBLAS_INLINE
1447	size_type one2 (size_type /* size1 /, size_type / size2 */) {
1448	return 1;
1449	}
1450
1451	static
1452	BOOST_UBLAS_INLINE
1453	size_type triangular_size (size_type size1, size_type size2) {
1454	size_type size = (std::max) (size1, size2);
1455	// Guard against size_type overflow - simplified
1456	BOOST_UBLAS_CHECK (size == 0 \|\| size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
1457	return ((size + 1) * size) / 2;
1458	}
1459	static
1460	BOOST_UBLAS_INLINE
1461	size_type lower_element (size_type i, size_type size1, size_type j, size_type size2) {
1462	BOOST_UBLAS_CHECK (i < size1, bad_index ());
1463	BOOST_UBLAS_CHECK (j < size2, bad_index ());
1464	BOOST_UBLAS_CHECK (i >= j, bad_index ());
1465	detail::ignore_unused_variable_warning(size1);
1466	detail::ignore_unused_variable_warning(size2);
1467	// FIXME size_type overflow
1468	// sigma_i (i + 1) = (i + 1) * i / 2
1469	// i = 0 1 2 3, sigma = 0 1 3 6
1470	return ((i + 1) * i) / 2 + j;
1471	}
1472	static
1473	BOOST_UBLAS_INLINE
1474	size_type upper_element (size_type i, size_type size1, size_type j, size_type size2) {
1475	BOOST_UBLAS_CHECK (i < size1, bad_index ());
1476	BOOST_UBLAS_CHECK (j < size2, bad_index ());
1477	BOOST_UBLAS_CHECK (i <= j, bad_index ());
1478	// FIXME size_type overflow
1479	// sigma_i (size - i) = size * i - i * (i - 1) / 2
1480	// i = 0 1 2 3, sigma = 0 4 7 9
1481	return (i * (2 * (std::max) (size1, size2) - i + 1)) / 2 + j - i;
1482	}
1483
1484	static
1485	BOOST_UBLAS_INLINE
1486	size_type element1 (size_type i, size_type size1, size_type /* j /, size_type / size2 */) {
1487	BOOST_UBLAS_CHECK (i < size1, bad_index ());
1488	detail::ignore_unused_variable_warning(size1);
1489	return i;
1490	}
1491	static
1492	BOOST_UBLAS_INLINE
1493	size_type element2 (size_type /* i /, size_type / size1 */, size_type j, size_type size2) {
1494	BOOST_UBLAS_CHECK (j < size2, bad_index ());
1495	detail::ignore_unused_variable_warning(size2);
1496	return j;
1497	}
1498	static
1499	BOOST_UBLAS_INLINE
1500	size_type address1 (size_type i, size_type size1, size_type /* j /, size_type / size2 */) {
1501	BOOST_UBLAS_CHECK (i <= size1, bad_index ());
1502	detail::ignore_unused_variable_warning(size1);
1503	return i;
1504	}
1505	static
1506	BOOST_UBLAS_INLINE
1507	size_type address2 (size_type /* i /, size_type / size1 */, size_type j, size_type size2) {
1508	BOOST_UBLAS_CHECK (j <= size2, bad_index ());
1509	detail::ignore_unused_variable_warning(size2);
1510	return j;
1511	}
1512	static
1513	BOOST_UBLAS_INLINE
1514	size_type index1 (size_type index1, size_type /* index2 */) {
1515	return index1;
1516	}
1517	static
1518	BOOST_UBLAS_INLINE
1519	size_type index2 (size_type /* index1 */, size_type index2) {
1520	return index2;
1521	}
1522	static
1523	BOOST_UBLAS_INLINE
1524	size_type size1 (size_type size1, size_type /* size2 */) {
1525	return size1;
1526	}
1527	static
1528	BOOST_UBLAS_INLINE
1529	size_type size2 (size_type /* size1 */, size_type size2) {
1530	return size2;
1531	}
1532
1533	// Iterating
1534	template<class I>
1535	static
1536	BOOST_UBLAS_INLINE
1537	void increment1 (I &it, size_type /* size1 */, size_type size2) {
1538	it += size2;
1539	}
1540	template<class I>
1541	static
1542	BOOST_UBLAS_INLINE
1543	void decrement1 (I &it, size_type /* size1 */, size_type size2) {
1544	it -= size2;
1545	}
1546	template<class I>
1547	static
1548	BOOST_UBLAS_INLINE
1549	void increment2 (I &it, size_type /* size1 /, size_type / size2 */) {
1550	++ it;
1551	}
1552	template<class I>
1553	static
1554	BOOST_UBLAS_INLINE
1555	void decrement2 (I &it, size_type /* size1 /, size_type / size2 */) {
1556	-- it;
1557	}
1558	};
1559
1560	// This functor computes the address translation
1561	// matrix [i] [j] -> storage [i + j * size1]
1562	template <class Z, class D>
1563	struct basic_column_major {
1564	typedef Z size_type;
1565	typedef D difference_type;
1566	typedef column_major_tag orientation_category;
1567
1568	static
1569	BOOST_UBLAS_INLINE
1570	size_type storage_size (size_type size1, size_type size2) {
1571	// Guard against size_type overflow
1572	BOOST_UBLAS_CHECK (size1 == 0 \|\| size2 <= (std::numeric_limits<size_type>::max) () / size1, bad_size ());
1573	return size1 * size2;
1574	}
1575
1576	// Indexing
1577	static
1578	BOOST_UBLAS_INLINE
1579	size_type element (size_type i, size_type size1, size_type j, size_type size2) {
1580	BOOST_UBLAS_CHECK (i < size1, bad_index ());
1581	BOOST_UBLAS_CHECK (j < size2, bad_index ());
1582	detail::ignore_unused_variable_warning(size2);
1583	// Guard against size_type overflow
1584	BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size1, bad_index ());
1585	return i + j * size1;
1586	}
1587	static
1588	BOOST_UBLAS_INLINE
1589	size_type address (size_type i, size_type size1, size_type j, size_type size2) {
1590	BOOST_UBLAS_CHECK (i <= size1, bad_index ());
1591	BOOST_UBLAS_CHECK (j <= size2, bad_index ());
1592	detail::ignore_unused_variable_warning(size2);
1593	// Guard against size_type overflow - address may be size1 past end of storage
1594	BOOST_UBLAS_CHECK (size1 == 0 \|\| j <= ((std::numeric_limits<size_type>::max) () - i) / size1, bad_index ());
1595	return i + j * size1;
1596	}
1597	static
1598	BOOST_UBLAS_INLINE
1599	difference_type distance1 (difference_type k, size_type /* size1 /, size_type / size2 */) {
1600	return k;
1601	}
1602	static
1603	BOOST_UBLAS_INLINE
1604	difference_type distance2 (difference_type k, size_type size1, size_type /* size2 */) {
1605	return size1 != 0 ? k / size1 : 0;
1606	}
1607	static
1608	BOOST_UBLAS_INLINE
1609	size_type index1 (difference_type k, size_type size1, size_type /* size2 */) {
1610	return size1 != 0 ? k % size1 : 0;
1611	}
1612	static
1613	BOOST_UBLAS_INLINE
1614	size_type index2 (difference_type k, size_type size1, size_type /* size2 */) {
1615	return size1 != 0 ? k / size1 : 0;
1616	}
1617	static
1618	BOOST_UBLAS_INLINE
1619	bool fast1 () {
1620	return true;
1621	}
1622	static
1623	BOOST_UBLAS_INLINE
1624	size_type one1 (size_type /* size1 /, size_type / size2 */) {
1625	return 1;
1626	}
1627	static
1628	BOOST_UBLAS_INLINE
1629	bool fast2 () {
1630	return false;
1631	}
1632	static
1633	BOOST_UBLAS_INLINE
1634	size_type one2 (size_type size1, size_type /* size2 */) {
1635	return size1;
1636	}
1637
1638	static
1639	BOOST_UBLAS_INLINE
1640	size_type triangular_size (size_type size1, size_type size2) {
1641	size_type size = (std::max) (size1, size2);
1642	// Guard against size_type overflow - simplified
1643	BOOST_UBLAS_CHECK (size == 0 \|\| size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
1644	return ((size + 1) * size) / 2;
1645	}
1646	static
1647	BOOST_UBLAS_INLINE
1648	size_type lower_element (size_type i, size_type size1, size_type j, size_type size2) {
1649	BOOST_UBLAS_CHECK (i < size1, bad_index ());
1650	BOOST_UBLAS_CHECK (j < size2, bad_index ());
1651	BOOST_UBLAS_CHECK (i >= j, bad_index ());
1652	// FIXME size_type overflow
1653	// sigma_j (size - j) = size * j - j * (j - 1) / 2
1654	// j = 0 1 2 3, sigma = 0 4 7 9
1655	return i - j + (j * (2 * (std::max) (size1, size2) - j + 1)) / 2;
1656	}
1657	static
1658	BOOST_UBLAS_INLINE
1659	size_type upper_element (size_type i, size_type size1, size_type j, size_type size2) {
1660	BOOST_UBLAS_CHECK (i < size1, bad_index ());
1661	BOOST_UBLAS_CHECK (j < size2, bad_index ());
1662	BOOST_UBLAS_CHECK (i <= j, bad_index ());
1663	// FIXME size_type overflow
1664	// sigma_j (j + 1) = (j + 1) * j / 2
1665	// j = 0 1 2 3, sigma = 0 1 3 6
1666	return i + ((j + 1) * j) / 2;
1667	}
1668
1669	static
1670	BOOST_UBLAS_INLINE
1671	size_type element1 (size_type /* i /, size_type / size1 */, size_type j, size_type size2) {
1672	BOOST_UBLAS_CHECK (j < size2, bad_index ());
1673	detail::ignore_unused_variable_warning(size2);
1674	return j;
1675	}
1676	static
1677	BOOST_UBLAS_INLINE
1678	size_type element2 (size_type i, size_type size1, size_type /* j /, size_type / size2 */) {
1679	BOOST_UBLAS_CHECK (i < size1, bad_index ());
1680	detail::ignore_unused_variable_warning(size1);
1681	return i;
1682	}
1683	static
1684	BOOST_UBLAS_INLINE
1685	size_type address1 (size_type /* i /, size_type / size1 */, size_type j, size_type size2) {
1686	BOOST_UBLAS_CHECK (j <= size2, bad_index ());
1687	detail::ignore_unused_variable_warning(size2);
1688	return j;
1689	}
1690	static
1691	BOOST_UBLAS_INLINE
1692	size_type address2 (size_type i, size_type size1, size_type /* j /, size_type / size2 */) {
1693	BOOST_UBLAS_CHECK (i <= size1, bad_index ());
1694	detail::ignore_unused_variable_warning(size1);
1695	return i;
1696	}
1697	static
1698	BOOST_UBLAS_INLINE
1699	size_type index1 (size_type /* index1 */, size_type index2) {
1700	return index2;
1701	}
1702	static
1703	BOOST_UBLAS_INLINE
1704	size_type index2 (size_type index1, size_type /* index2 */) {
1705	return index1;
1706	}
1707	static
1708	BOOST_UBLAS_INLINE
1709	size_type size1 (size_type /* size1 */, size_type size2) {
1710	return size2;
1711	}
1712	static
1713	BOOST_UBLAS_INLINE
1714	size_type size2 (size_type size1, size_type /* size2 */) {
1715	return size1;
1716	}
1717
1718	// Iterating
1719	template<class I>
1720	static
1721	BOOST_UBLAS_INLINE
1722	void increment1 (I &it, size_type /* size1 /, size_type / size2 */) {
1723	++ it;
1724	}
1725	template<class I>
1726	static
1727	BOOST_UBLAS_INLINE
1728	void decrement1 (I &it, size_type /* size1 /, size_type / size2 */) {
1729	-- it;
1730	}
1731	template<class I>
1732	static
1733	BOOST_UBLAS_INLINE
1734	void increment2 (I &it, size_type size1, size_type /* size2 */) {
1735	it += size1;
1736	}
1737	template<class I>
1738	static
1739	BOOST_UBLAS_INLINE
1740	void decrement2 (I &it, size_type size1, size_type /* size2 */) {
1741	it -= size1;
1742	}
1743	};
1744
1745
1746	template <class Z>
1747	struct basic_full {
1748	typedef Z size_type;
1749
1750	template<class L>
1751	static
1752	BOOST_UBLAS_INLINE
1753	size_type packed_size (size_type size1, size_type size2) {
1754	return L::storage_size (size1, size2);
1755	}
1756
1757	static
1758	BOOST_UBLAS_INLINE
1759	bool zero (size_type /* i /, size_type / j */) {
1760	return false;
1761	}
1762	static
1763	BOOST_UBLAS_INLINE
1764	bool one (size_type /* i /, size_type / j */) {
1765	return false;
1766	}
1767	static
1768	BOOST_UBLAS_INLINE
1769	bool other (size_type /* i /, size_type / j */) {
1770	return true;
1771	}
1772	};
1773
1774	template <class Z>
1775	struct basic_lower {
1776	typedef Z size_type;
1777
1778	template<class L>
1779	static
1780	BOOST_UBLAS_INLINE
1781	size_type packed_size (L, size_type size1, size_type size2) {
1782	return L::triangular_size (size1, size2);
1783	}
1784
1785	static
1786	BOOST_UBLAS_INLINE
1787	bool zero (size_type i, size_type j) {
1788	return j > i;
1789	}
1790	static
1791	BOOST_UBLAS_INLINE
1792	bool one (size_type /* i /, size_type / j */) {
1793	return false;
1794	}
1795	static
1796	BOOST_UBLAS_INLINE
1797	bool other (size_type i, size_type j) {
1798	return j <= i;
1799	}
1800	template<class L>
1801	static
1802	BOOST_UBLAS_INLINE
1803	size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1804	return L::lower_element (i, size1, j, size2);
1805	}
1806
1807	static
1808	BOOST_UBLAS_INLINE
1809	size_type restrict1 (size_type i, size_type j) {
1810	return (std::max) (i, j);
1811	}
1812	static
1813	BOOST_UBLAS_INLINE
1814	size_type restrict2 (size_type i, size_type j) {
1815	return (std::min) (i + 1, j);
1816	}
1817	static
1818	BOOST_UBLAS_INLINE
1819	size_type mutable_restrict1 (size_type i, size_type j) {
1820	return (std::max) (i, j);
1821	}
1822	static
1823	BOOST_UBLAS_INLINE
1824	size_type mutable_restrict2 (size_type i, size_type j) {
1825	return (std::min) (i + 1, j);
1826	}
1827	};
1828	template <class Z>
1829	struct basic_upper {
1830	typedef Z size_type;
1831
1832	template<class L>
1833	static
1834	BOOST_UBLAS_INLINE
1835	size_type packed_size (L, size_type size1, size_type size2) {
1836	return L::triangular_size (size1, size2);
1837	}
1838
1839	static
1840	BOOST_UBLAS_INLINE
1841	bool zero (size_type i, size_type j) {
1842	return j < i;
1843	}
1844	static
1845	BOOST_UBLAS_INLINE
1846	bool one (size_type /* i /, size_type / j */) {
1847	return false;
1848	}
1849	static
1850	BOOST_UBLAS_INLINE
1851	bool other (size_type i, size_type j) {
1852	return j >= i;
1853	}
1854	template<class L>
1855	static
1856	BOOST_UBLAS_INLINE
1857	size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1858	return L::upper_element (i, size1, j, size2);
1859	}
1860
1861	static
1862	BOOST_UBLAS_INLINE
1863	size_type restrict1 (size_type i, size_type j) {
1864	return (std::min) (i, j + 1);
1865	}
1866	static
1867	BOOST_UBLAS_INLINE
1868	size_type restrict2 (size_type i, size_type j) {
1869	return (std::max) (i, j);
1870	}
1871	static
1872	BOOST_UBLAS_INLINE
1873	size_type mutable_restrict1 (size_type i, size_type j) {
1874	return (std::min) (i, j + 1);
1875	}
1876	static
1877	BOOST_UBLAS_INLINE
1878	size_type mutable_restrict2 (size_type i, size_type j) {
1879	return (std::max) (i, j);
1880	}
1881	};
1882	template <class Z>
1883	struct basic_unit_lower : public basic_lower<Z> {
1884	typedef Z size_type;
1885
1886	template<class L>
1887	static
1888	BOOST_UBLAS_INLINE
1889	size_type packed_size (L, size_type size1, size_type size2) {
1890	// Zero size strict triangles are bad at this point
1891	BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1892	return L::triangular_size (size1 - 1, size2 - 1);
1893	}
1894
1895	static
1896	BOOST_UBLAS_INLINE
1897	bool one (size_type i, size_type j) {
1898	return j == i;
1899	}
1900	static
1901	BOOST_UBLAS_INLINE
1902	bool other (size_type i, size_type j) {
1903	return j < i;
1904	}
1905	template<class L>
1906	static
1907	BOOST_UBLAS_INLINE
1908	size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1909	// Zero size strict triangles are bad at this point
1910	BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1911	return L::lower_element (i, size1 - 1, j, size2 - 1);
1912	}
1913
1914	static
1915	BOOST_UBLAS_INLINE
1916	size_type mutable_restrict1 (size_type i, size_type j) {
1917	return (std::max) (i, j);
1918	}
1919	static
1920	BOOST_UBLAS_INLINE
1921	size_type mutable_restrict2 (size_type i, size_type j) {
1922	return (std::min) (i, j);
1923	}
1924	};
1925	template <class Z>
1926	struct basic_unit_upper : public basic_upper<Z> {
1927	typedef Z size_type;
1928
1929	template<class L>
1930	static
1931	BOOST_UBLAS_INLINE
1932	size_type packed_size (L, size_type size1, size_type size2) {
1933	// Zero size strict triangles are bad at this point
1934	BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1935	return L::triangular_size (size1 - 1, size2 - 1);
1936	}
1937
1938	static
1939	BOOST_UBLAS_INLINE
1940	bool one (size_type i, size_type j) {
1941	return j == i;
1942	}
1943	static
1944	BOOST_UBLAS_INLINE
1945	bool other (size_type i, size_type j) {
1946	return j > i;
1947	}
1948	template<class L>
1949	static
1950	BOOST_UBLAS_INLINE
1951	size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1952	// Zero size strict triangles are bad at this point
1953	BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1954	return L::upper_element (i, size1 - 1, j, size2 - 1);
1955	}
1956
1957	static
1958	BOOST_UBLAS_INLINE
1959	size_type mutable_restrict1 (size_type i, size_type j) {
1960	return (std::min) (i, j);
1961	}
1962	static
1963	BOOST_UBLAS_INLINE
1964	size_type mutable_restrict2 (size_type i, size_type j) {
1965	return (std::max) (i, j);
1966	}
1967	};
1968	template <class Z>
1969	struct basic_strict_lower : public basic_lower<Z> {
1970	typedef Z size_type;
1971
1972	template<class L>
1973	static
1974	BOOST_UBLAS_INLINE
1975	size_type packed_size (L, size_type size1, size_type size2) {
1976	// Zero size strict triangles are bad at this point
1977	BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1978	return L::triangular_size (size1 - 1, size2 - 1);
1979	}
1980
1981	static
1982	BOOST_UBLAS_INLINE
1983	bool zero (size_type i, size_type j) {
1984	return j >= i;
1985	}
1986	static
1987	BOOST_UBLAS_INLINE
1988	bool one (size_type /i/, size_type /j/) {
1989	return false;
1990	}
1991	static
1992	BOOST_UBLAS_INLINE
1993	bool other (size_type i, size_type j) {
1994	return j < i;
1995	}
1996	template<class L>
1997	static
1998	BOOST_UBLAS_INLINE
1999	size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
2000	// Zero size strict triangles are bad at this point
2001	BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
2002	return L::lower_element (i, size1 - 1, j, size2 - 1);
2003	}
2004
2005	static
2006	BOOST_UBLAS_INLINE
2007	size_type restrict1 (size_type i, size_type j) {
2008	return (std::max) (i, j);
2009	}
2010	static
2011	BOOST_UBLAS_INLINE
2012	size_type restrict2 (size_type i, size_type j) {
2013	return (std::min) (i, j);
2014	}
2015	static
2016	BOOST_UBLAS_INLINE
2017	size_type mutable_restrict1 (size_type i, size_type j) {
2018	return (std::max) (i, j);
2019	}
2020	static
2021	BOOST_UBLAS_INLINE
2022	size_type mutable_restrict2 (size_type i, size_type j) {
2023	return (std::min) (i, j);
2024	}
2025	};
2026	template <class Z>
2027	struct basic_strict_upper : public basic_upper<Z> {
2028	typedef Z size_type;
2029
2030	template<class L>
2031	static
2032	BOOST_UBLAS_INLINE
2033	size_type packed_size (L, size_type size1, size_type size2) {
2034	// Zero size strict triangles are bad at this point
2035	BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
2036	return L::triangular_size (size1 - 1, size2 - 1);
2037	}
2038
2039	static
2040	BOOST_UBLAS_INLINE
2041	bool zero (size_type i, size_type j) {
2042	return j <= i;
2043	}
2044	static
2045	BOOST_UBLAS_INLINE
2046	bool one (size_type /i/, size_type /j/) {
2047	return false;
2048	}
2049	static
2050	BOOST_UBLAS_INLINE
2051	bool other (size_type i, size_type j) {
2052	return j > i;
2053	}
2054	template<class L>
2055	static
2056	BOOST_UBLAS_INLINE
2057	size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
2058	// Zero size strict triangles are bad at this point
2059	BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
2060	return L::upper_element (i, size1 - 1, j, size2 - 1);
2061	}
2062
2063	static
2064	BOOST_UBLAS_INLINE
2065	size_type restrict1 (size_type i, size_type j) {
2066	return (std::min) (i, j);
2067	}
2068	static
2069	BOOST_UBLAS_INLINE
2070	size_type restrict2 (size_type i, size_type j) {
2071	return (std::max) (i, j);
2072	}
2073	static
2074	BOOST_UBLAS_INLINE
2075	size_type mutable_restrict1 (size_type i, size_type j) {
2076	return (std::min) (i, j);
2077	}
2078	static
2079	BOOST_UBLAS_INLINE
2080	size_type mutable_restrict2 (size_type i, size_type j) {
2081	return (std::max) (i, j);
2082	}
2083	};
2084
2085	}}}
2086
2087	#endif

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

Context Navigation

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

Download in other formats: