Symaptic: os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/functional.hpp@bde4ae8d615e (annotated)

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

author	sl@SLION-WIN7.fritz.box
	Fri, 15 Jun 2012 03:10:57 +0200
changeset 0	bde4ae8d615e
permissions	-rw-r--r--