Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

functional.hpp @ 857

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

Rev	Line
[857]	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: