Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

operation.hpp @ 857

Revision 857, 33.5 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_OPERATION_
	18	#define _BOOST_UBLAS_OPERATION_
	19
	20	#include <boost/numeric/ublas/matrix_proxy.hpp>
	21
	22	/** \file operation.hpp
	23	* \brief This file contains some specialized products.
	24	*/
	25
	26	// axpy-based products
	27	// Alexei Novakov had a lot of ideas to improve these. Thanks.
	28	// Hendrik Kueck proposed some new kernel. Thanks again.
	29
	30	namespace boost { namespace numeric { namespace ublas {
	31
	32	template<class V, class T1, class IA1, class TA1, class E2>
	33	BOOST_UBLAS_INLINE
	34	V &
	35	axpy_prod (const compressed_matrix<T1, row_major, 0, IA1, TA1> &e1,
	36	const vector_expression<E2> &e2,
	37	V &v, row_major_tag) {
	38	typedef typename V::size_type size_type;
	39	typedef typename V::value_type value_type;
	40
	41	for (size_type i = 0; i < e1.filled1 () -1; ++ i) {
	42	size_type begin = e1.index1_data () [i];
	43	size_type end = e1.index1_data () [i + 1];
	44	value_type t (v (i));
	45	for (size_type j = begin; j < end; ++ j)
	46	t += e1.value_data () [j] * e2 () (e1.index2_data () [j]);
	47	v (i) = t;
	48	}
	49	return v;
	50	}
	51
	52	template<class V, class T1, class IA1, class TA1, class E2>
	53	BOOST_UBLAS_INLINE
	54	V &
	55	axpy_prod (const compressed_matrix<T1, column_major, 0, IA1, TA1> &e1,
	56	const vector_expression<E2> &e2,
	57	V &v, column_major_tag) {
	58	typedef typename V::size_type size_type;
	59
	60	for (size_type j = 0; j < e1.filled1 () -1; ++ j) {
	61	size_type begin = e1.index1_data () [j];
	62	size_type end = e1.index1_data () [j + 1];
	63	for (size_type i = begin; i < end; ++ i)
	64	v (e1.index2_data () [i]) += e1.value_data () [i] * e2 () (j);
	65	}
	66	return v;
	67	}
	68
	69	// Dispatcher
	70	template<class V, class T1, class L1, class IA1, class TA1, class E2>
	71	BOOST_UBLAS_INLINE
	72	V &
	73	axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
	74	const vector_expression<E2> &e2,
	75	V &v, bool init = true) {
	76	typedef typename V::value_type value_type;
	77	typedef typename L1::orientation_category orientation_category;
	78
	79	if (init)
	80	v.assign (zero_vector<value_type> (e1.size1 ()));
	81	#if BOOST_UBLAS_TYPE_CHECK
	82	vector<value_type> cv (v);
	83	typedef typename type_traits<value_type>::real_type real_type;
	84	real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
	85	indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
	86	#endif
	87	axpy_prod (e1, e2, v, orientation_category ());
	88	#if BOOST_UBLAS_TYPE_CHECK
	89	BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
	90	#endif
	91	return v;
	92	}
	93	template<class V, class T1, class L1, class IA1, class TA1, class E2>
	94	BOOST_UBLAS_INLINE
	95	V
	96	axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
	97	const vector_expression<E2> &e2) {
	98	typedef V vector_type;
	99
	100	vector_type v (e1.size1 ());
	101	return axpy_prod (e1, e2, v, true);
	102	}
	103
	104	template<class V, class T1, class L1, class IA1, class TA1, class E2>
	105	BOOST_UBLAS_INLINE
	106	V &
	107	axpy_prod (const coordinate_matrix<T1, L1, 0, IA1, TA1> &e1,
	108	const vector_expression<E2> &e2,
	109	V &v, bool init = true) {
	110	typedef typename V::size_type size_type;
	111	typedef typename V::value_type value_type;
	112	typedef L1 layout_type;
	113
	114	size_type size1 = e1.size1();
	115	size_type size2 = e1.size2();
	116
	117	if (init) {
	118	noalias(v) = zero_vector<value_type>(size1);
	119	}
	120
	121	for (size_type i = 0; i < e1.nnz(); ++i) {
	122	size_type row_index = layout_type::element1( e1.index1_data () [i], size1, e1.index2_data () [i], size2 );
	123	size_type col_index = layout_type::element2( e1.index1_data () [i], size1, e1.index2_data () [i], size2 );
	124	v( row_index ) += e1.value_data () [i] * e2 () (col_index);
	125	}
	126	return v;
	127	}
	128
	129	template<class V, class E1, class E2>
	130	BOOST_UBLAS_INLINE
	131	V &
	132	axpy_prod (const matrix_expression<E1> &e1,
	133	const vector_expression<E2> &e2,
	134	V &v, packed_random_access_iterator_tag, row_major_tag) {
	135	typedef const E1 expression1_type;
	136	typedef const E2 expression2_type;
	137	typedef typename V::size_type size_type;
	138
	139	typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
	140	typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
	141	while (it1 != it1_end) {
	142	size_type index1 (it1.index1 ());
	143	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
	144	typename expression1_type::const_iterator2 it2 (it1.begin ());
	145	typename expression1_type::const_iterator2 it2_end (it1.end ());
	146	#else
	147	typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
	148	typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
	149	#endif
	150	while (it2 != it2_end) {
	151	v (index1) += it2 e2 () (it2.index2 ());
	152	++ it2;
	153	}
	154	++ it1;
	155	}
	156	return v;
	157	}
	158
	159	template<class V, class E1, class E2>
	160	BOOST_UBLAS_INLINE
	161	V &
	162	axpy_prod (const matrix_expression<E1> &e1,
	163	const vector_expression<E2> &e2,
	164	V &v, packed_random_access_iterator_tag, column_major_tag) {
	165	typedef const E1 expression1_type;
	166	typedef const E2 expression2_type;
	167	typedef typename V::size_type size_type;
	168
	169	typename expression1_type::const_iterator2 it2 (e1 ().begin2 ());
	170	typename expression1_type::const_iterator2 it2_end (e1 ().end2 ());
	171	while (it2 != it2_end) {
	172	size_type index2 (it2.index2 ());
	173	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
	174	typename expression1_type::const_iterator1 it1 (it2.begin ());
	175	typename expression1_type::const_iterator1 it1_end (it2.end ());
	176	#else
	177	typename expression1_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
	178	typename expression1_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
	179	#endif
	180	while (it1 != it1_end) {
	181	v (it1.index1 ()) += it1 e2 () (index2);
	182	++ it1;
	183	}
	184	++ it2;
	185	}
	186	return v;
	187	}
	188
	189	template<class V, class E1, class E2>
	190	BOOST_UBLAS_INLINE
	191	V &
	192	axpy_prod (const matrix_expression<E1> &e1,
	193	const vector_expression<E2> &e2,
	194	V &v, sparse_bidirectional_iterator_tag) {
	195	typedef const E1 expression1_type;
	196	typedef const E2 expression2_type;
	197	typedef typename V::size_type size_type;
	198
	199	typename expression2_type::const_iterator it (e2 ().begin ());
	200	typename expression2_type::const_iterator it_end (e2 ().end ());
	201	while (it != it_end) {
	202	v.plus_assign (column (e1 (), it.index ()) * *it);
	203	++ it;
	204	}
	205	return v;
	206	}
	207
	208	// Dispatcher
	209	template<class V, class E1, class E2>
	210	BOOST_UBLAS_INLINE
	211	V &
	212	axpy_prod (const matrix_expression<E1> &e1,
	213	const vector_expression<E2> &e2,
	214	V &v, packed_random_access_iterator_tag) {
	215	typedef typename E1::orientation_category orientation_category;
	216	return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
	217	}
	218
	219
	220	/** \brief computes <tt>v += A x</tt> or <tt>v = A x</tt> in an
	221	optimized fashion.
	222
	223	\param e1 the matrix expression \c A
	224	\param e2 the vector expression \c x
	225	\param v the result vector \c v
	226	\param init a boolean parameter
	227
	228	<tt>axpy_prod(A, x, v, init)</tt> implements the well known
	229	axpy-product. Setting \a init to \c true is equivalent to call
	230	<tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
	231	defaults to \c true, but this may change in the future.
	232
	233	Up to now there are some specialisation for compressed
	234	matrices that give a large speed up compared to prod.
	235
	236	\ingroup blas2
	237
	238	\internal
	239
	240	template parameters:
	241	\param V type of the result vector \c v
	242	\param E1 type of a matrix expression \c A
	243	\param E2 type of a vector expression \c x
	244	*/
	245	template<class V, class E1, class E2>
	246	BOOST_UBLAS_INLINE
	247	V &
	248	axpy_prod (const matrix_expression<E1> &e1,
	249	const vector_expression<E2> &e2,
	250	V &v, bool init = true) {
	251	typedef typename V::value_type value_type;
	252	typedef typename E2::const_iterator::iterator_category iterator_category;
	253
	254	if (init)
	255	v.assign (zero_vector<value_type> (e1 ().size1 ()));
	256	#if BOOST_UBLAS_TYPE_CHECK
	257	vector<value_type> cv (v);
	258	typedef typename type_traits<value_type>::real_type real_type;
	259	real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
	260	indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
	261	#endif
	262	axpy_prod (e1, e2, v, iterator_category ());
	263	#if BOOST_UBLAS_TYPE_CHECK
	264	BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
	265	#endif
	266	return v;
	267	}
	268	template<class V, class E1, class E2>
	269	BOOST_UBLAS_INLINE
	270	V
	271	axpy_prod (const matrix_expression<E1> &e1,
	272	const vector_expression<E2> &e2) {
	273	typedef V vector_type;
	274
	275	vector_type v (e1 ().size1 ());
	276	return axpy_prod (e1, e2, v, true);
	277	}
	278
	279	template<class V, class E1, class T2, class IA2, class TA2>
	280	BOOST_UBLAS_INLINE
	281	V &
	282	axpy_prod (const vector_expression<E1> &e1,
	283	const compressed_matrix<T2, column_major, 0, IA2, TA2> &e2,
	284	V &v, column_major_tag) {
	285	typedef typename V::size_type size_type;
	286	typedef typename V::value_type value_type;
	287
	288	for (size_type j = 0; j < e2.filled1 () -1; ++ j) {
	289	size_type begin = e2.index1_data () [j];
	290	size_type end = e2.index1_data () [j + 1];
	291	value_type t (v (j));
	292	for (size_type i = begin; i < end; ++ i)
	293	t += e2.value_data () [i] * e1 () (e2.index2_data () [i]);
	294	v (j) = t;
	295	}
	296	return v;
	297	}
	298
	299	template<class V, class E1, class T2, class IA2, class TA2>
	300	BOOST_UBLAS_INLINE
	301	V &
	302	axpy_prod (const vector_expression<E1> &e1,
	303	const compressed_matrix<T2, row_major, 0, IA2, TA2> &e2,
	304	V &v, row_major_tag) {
	305	typedef typename V::size_type size_type;
	306
	307	for (size_type i = 0; i < e2.filled1 () -1; ++ i) {
	308	size_type begin = e2.index1_data () [i];
	309	size_type end = e2.index1_data () [i + 1];
	310	for (size_type j = begin; j < end; ++ j)
	311	v (e2.index2_data () [j]) += e2.value_data () [j] * e1 () (i);
	312	}
	313	return v;
	314	}
	315
	316	// Dispatcher
	317	template<class V, class E1, class T2, class L2, class IA2, class TA2>
	318	BOOST_UBLAS_INLINE
	319	V &
	320	axpy_prod (const vector_expression<E1> &e1,
	321	const compressed_matrix<T2, L2, 0, IA2, TA2> &e2,
	322	V &v, bool init = true) {
	323	typedef typename V::value_type value_type;
	324	typedef typename L2::orientation_category orientation_category;
	325
	326	if (init)
	327	v.assign (zero_vector<value_type> (e2.size2 ()));
	328	#if BOOST_UBLAS_TYPE_CHECK
	329	vector<value_type> cv (v);
	330	typedef typename type_traits<value_type>::real_type real_type;
	331	real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
	332	indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
	333	#endif
	334	axpy_prod (e1, e2, v, orientation_category ());
	335	#if BOOST_UBLAS_TYPE_CHECK
	336	BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
	337	#endif
	338	return v;
	339	}
	340	template<class V, class E1, class T2, class L2, class IA2, class TA2>
	341	BOOST_UBLAS_INLINE
	342	V
	343	axpy_prod (const vector_expression<E1> &e1,
	344	const compressed_matrix<T2, L2, 0, IA2, TA2> &e2) {
	345	typedef V vector_type;
	346
	347	vector_type v (e2.size2 ());
	348	return axpy_prod (e1, e2, v, true);
	349	}
	350
	351	template<class V, class E1, class E2>
	352	BOOST_UBLAS_INLINE
	353	V &
	354	axpy_prod (const vector_expression<E1> &e1,
	355	const matrix_expression<E2> &e2,
	356	V &v, packed_random_access_iterator_tag, column_major_tag) {
	357	typedef const E1 expression1_type;
	358	typedef const E2 expression2_type;
	359	typedef typename V::size_type size_type;
	360
	361	typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
	362	typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
	363	while (it2 != it2_end) {
	364	size_type index2 (it2.index2 ());
	365	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
	366	typename expression2_type::const_iterator1 it1 (it2.begin ());
	367	typename expression2_type::const_iterator1 it1_end (it2.end ());
	368	#else
	369	typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
	370	typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
	371	#endif
	372	while (it1 != it1_end) {
	373	v (index2) += it1 e1 () (it1.index1 ());
	374	++ it1;
	375	}
	376	++ it2;
	377	}
	378	return v;
	379	}
	380
	381	template<class V, class E1, class E2>
	382	BOOST_UBLAS_INLINE
	383	V &
	384	axpy_prod (const vector_expression<E1> &e1,
	385	const matrix_expression<E2> &e2,
	386	V &v, packed_random_access_iterator_tag, row_major_tag) {
	387	typedef const E1 expression1_type;
	388	typedef const E2 expression2_type;
	389	typedef typename V::size_type size_type;
	390
	391	typename expression2_type::const_iterator1 it1 (e2 ().begin1 ());
	392	typename expression2_type::const_iterator1 it1_end (e2 ().end1 ());
	393	while (it1 != it1_end) {
	394	size_type index1 (it1.index1 ());
	395	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
	396	typename expression2_type::const_iterator2 it2 (it1.begin ());
	397	typename expression2_type::const_iterator2 it2_end (it1.end ());
	398	#else
	399	typename expression2_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
	400	typename expression2_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
	401	#endif
	402	while (it2 != it2_end) {
	403	v (it2.index2 ()) += it2 e1 () (index1);
	404	++ it2;
	405	}
	406	++ it1;
	407	}
	408	return v;
	409	}
	410
	411	template<class V, class E1, class E2>
	412	BOOST_UBLAS_INLINE
	413	V &
	414	axpy_prod (const vector_expression<E1> &e1,
	415	const matrix_expression<E2> &e2,
	416	V &v, sparse_bidirectional_iterator_tag) {
	417	typedef const E1 expression1_type;
	418	typedef const E2 expression2_type;
	419	typedef typename V::size_type size_type;
	420
	421	typename expression1_type::const_iterator it (e1 ().begin ());
	422	typename expression1_type::const_iterator it_end (e1 ().end ());
	423	while (it != it_end) {
	424	v.plus_assign (it row (e2 (), it.index ()));
	425	++ it;
	426	}
	427	return v;
	428	}
	429
	430	// Dispatcher
	431	template<class V, class E1, class E2>
	432	BOOST_UBLAS_INLINE
	433	V &
	434	axpy_prod (const vector_expression<E1> &e1,
	435	const matrix_expression<E2> &e2,
	436	V &v, packed_random_access_iterator_tag) {
	437	typedef typename E2::orientation_category orientation_category;
	438	return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
	439	}
	440
	441
	442	/** \brief computes <tt>v += A<sup>T</sup> x</tt> or <tt>v = A<sup>T</sup> x</tt> in an
	443	optimized fashion.
	444
	445	\param e1 the vector expression \c x
	446	\param e2 the matrix expression \c A
	447	\param v the result vector \c v
	448	\param init a boolean parameter
	449
	450	<tt>axpy_prod(x, A, v, init)</tt> implements the well known
	451	axpy-product. Setting \a init to \c true is equivalent to call
	452	<tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
	453	defaults to \c true, but this may change in the future.
	454
	455	Up to now there are some specialisation for compressed
	456	matrices that give a large speed up compared to prod.
	457
	458	\ingroup blas2
	459
	460	\internal
	461
	462	template parameters:
	463	\param V type of the result vector \c v
	464	\param E1 type of a vector expression \c x
	465	\param E2 type of a matrix expression \c A
	466	*/
	467	template<class V, class E1, class E2>
	468	BOOST_UBLAS_INLINE
	469	V &
	470	axpy_prod (const vector_expression<E1> &e1,
	471	const matrix_expression<E2> &e2,
	472	V &v, bool init = true) {
	473	typedef typename V::value_type value_type;
	474	typedef typename E1::const_iterator::iterator_category iterator_category;
	475
	476	if (init)
	477	v.assign (zero_vector<value_type> (e2 ().size2 ()));
	478	#if BOOST_UBLAS_TYPE_CHECK
	479	vector<value_type> cv (v);
	480	typedef typename type_traits<value_type>::real_type real_type;
	481	real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
	482	indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
	483	#endif
	484	axpy_prod (e1, e2, v, iterator_category ());
	485	#if BOOST_UBLAS_TYPE_CHECK
	486	BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
	487	#endif
	488	return v;
	489	}
	490	template<class V, class E1, class E2>
	491	BOOST_UBLAS_INLINE
	492	V
	493	axpy_prod (const vector_expression<E1> &e1,
	494	const matrix_expression<E2> &e2) {
	495	typedef V vector_type;
	496
	497	vector_type v (e2 ().size2 ());
	498	return axpy_prod (e1, e2, v, true);
	499	}
	500
	501	template<class M, class E1, class E2, class TRI>
	502	BOOST_UBLAS_INLINE
	503	M &
	504	axpy_prod (const matrix_expression<E1> &e1,
	505	const matrix_expression<E2> &e2,
	506	M &m, TRI,
	507	dense_proxy_tag, row_major_tag) {
	508	typedef M matrix_type;
	509	typedef const E1 expression1_type;
	510	typedef const E2 expression2_type;
	511	typedef typename M::size_type size_type;
	512	typedef typename M::value_type value_type;
	513
	514	#if BOOST_UBLAS_TYPE_CHECK
	515	matrix<value_type, row_major> cm (m);
	516	typedef typename type_traits<value_type>::real_type real_type;
	517	real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
	518	indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
	519	#endif
	520	size_type size1 (e1 ().size1 ());
	521	size_type size2 (e1 ().size2 ());
	522	for (size_type i = 0; i < size1; ++ i)
	523	for (size_type j = 0; j < size2; ++ j)
	524	row (m, i).plus_assign (e1 () (i, j) * row (e2 (), j));
	525	#if BOOST_UBLAS_TYPE_CHECK
	526	BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
	527	#endif
	528	return m;
	529	}
	530	template<class M, class E1, class E2, class TRI>
	531	BOOST_UBLAS_INLINE
	532	M &
	533	axpy_prod (const matrix_expression<E1> &e1,
	534	const matrix_expression<E2> &e2,
	535	M &m, TRI,
	536	sparse_proxy_tag, row_major_tag) {
	537	typedef M matrix_type;
	538	typedef TRI triangular_restriction;
	539	typedef const E1 expression1_type;
	540	typedef const E2 expression2_type;
	541	typedef typename M::size_type size_type;
	542	typedef typename M::value_type value_type;
	543
	544	#if BOOST_UBLAS_TYPE_CHECK
	545	matrix<value_type, row_major> cm (m);
	546	typedef typename type_traits<value_type>::real_type real_type;
	547	real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
	548	indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
	549	#endif
	550	typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
	551	typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
	552	while (it1 != it1_end) {
	553	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
	554	typename expression1_type::const_iterator2 it2 (it1.begin ());
	555	typename expression1_type::const_iterator2 it2_end (it1.end ());
	556	#else
	557	typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
	558	typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
	559	#endif
	560	while (it2 != it2_end) {
	561	// row (m, it1.index1 ()).plus_assign (it2 row (e2 (), it2.index2 ()));
	562	matrix_row<expression2_type> mr (e2 (), it2.index2 ());
	563	typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
	564	typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
	565	while (itr != itr_end) {
	566	if (triangular_restriction::other (it1.index1 (), itr.index ()))
	567	m (it1.index1 (), itr.index ()) += it2 *itr;
	568	++ itr;
	569	}
	570	++ it2;
	571	}
	572	++ it1;
	573	}
	574	#if BOOST_UBLAS_TYPE_CHECK
	575	BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
	576	#endif
	577	return m;
	578	}
	579
	580	template<class M, class E1, class E2, class TRI>
	581	BOOST_UBLAS_INLINE
	582	M &
	583	axpy_prod (const matrix_expression<E1> &e1,
	584	const matrix_expression<E2> &e2,
	585	M &m, TRI,
	586	dense_proxy_tag, column_major_tag) {
	587	typedef M matrix_type;
	588	typedef const E1 expression1_type;
	589	typedef const E2 expression2_type;
	590	typedef typename M::size_type size_type;
	591	typedef typename M::value_type value_type;
	592
	593	#if BOOST_UBLAS_TYPE_CHECK
	594	matrix<value_type, column_major> cm (m);
	595	typedef typename type_traits<value_type>::real_type real_type;
	596	real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
	597	indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
	598	#endif
	599	size_type size1 (e2 ().size1 ());
	600	size_type size2 (e2 ().size2 ());
	601	for (size_type j = 0; j < size2; ++ j)
	602	for (size_type i = 0; i < size1; ++ i)
	603	column (m, j).plus_assign (e2 () (i, j) * column (e1 (), i));
	604	#if BOOST_UBLAS_TYPE_CHECK
	605	BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
	606	#endif
	607	return m;
	608	}
	609	template<class M, class E1, class E2, class TRI>
	610	BOOST_UBLAS_INLINE
	611	M &
	612	axpy_prod (const matrix_expression<E1> &e1,
	613	const matrix_expression<E2> &e2,
	614	M &m, TRI,
	615	sparse_proxy_tag, column_major_tag) {
	616	typedef M matrix_type;
	617	typedef TRI triangular_restriction;
	618	typedef const E1 expression1_type;
	619	typedef const E2 expression2_type;
	620	typedef typename M::size_type size_type;
	621	typedef typename M::value_type value_type;
	622
	623	#if BOOST_UBLAS_TYPE_CHECK
	624	matrix<value_type, column_major> cm (m);
	625	typedef typename type_traits<value_type>::real_type real_type;
	626	real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
	627	indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
	628	#endif
	629	typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
	630	typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
	631	while (it2 != it2_end) {
	632	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
	633	typename expression2_type::const_iterator1 it1 (it2.begin ());
	634	typename expression2_type::const_iterator1 it1_end (it2.end ());
	635	#else
	636	typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
	637	typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
	638	#endif
	639	while (it1 != it1_end) {
	640	// column (m, it2.index2 ()).plus_assign (it1 column (e1 (), it1.index1 ()));
	641	matrix_column<expression1_type> mc (e1 (), it1.index1 ());
	642	typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
	643	typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
	644	while (itc != itc_end) {
	645	if (triangular_restriction::functor_type ().other (itc.index (), it2.index2 ()))
	646	m (itc.index (), it2.index2 ()) += it1 *itc;
	647	++ itc;
	648	}
	649	++ it1;
	650	}
	651	++ it2;
	652	}
	653	#if BOOST_UBLAS_TYPE_CHECK
	654	BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
	655	#endif
	656	return m;
	657	}
	658
	659	// Dispatcher
	660	template<class M, class E1, class E2, class TRI>
	661	BOOST_UBLAS_INLINE
	662	M &
	663	axpy_prod (const matrix_expression<E1> &e1,
	664	const matrix_expression<E2> &e2,
	665	M &m, TRI, bool init = true) {
	666	typedef typename M::value_type value_type;
	667	typedef typename M::storage_category storage_category;
	668	typedef typename M::orientation_category orientation_category;
	669	typedef TRI triangular_restriction;
	670
	671	if (init)
	672	m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
	673	return axpy_prod (e1, e2, m, triangular_restriction (), storage_category (), orientation_category ());
	674	}
	675	template<class M, class E1, class E2, class TRI>
	676	BOOST_UBLAS_INLINE
	677	M
	678	axpy_prod (const matrix_expression<E1> &e1,
	679	const matrix_expression<E2> &e2,
	680	TRI) {
	681	typedef M matrix_type;
	682	typedef TRI triangular_restriction;
	683
	684	matrix_type m (e1 ().size1 (), e2 ().size2 ());
	685	return axpy_prod (e1, e2, m, triangular_restriction (), true);
	686	}
	687
	688	/** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
	689	optimized fashion.
	690
	691	\param e1 the matrix expression \c A
	692	\param e2 the matrix expression \c X
	693	\param m the result matrix \c M
	694	\param init a boolean parameter
	695
	696	<tt>axpy_prod(A, X, M, init)</tt> implements the well known
	697	axpy-product. Setting \a init to \c true is equivalent to call
	698	<tt>M.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
	699	defaults to \c true, but this may change in the future.
	700
	701	Up to now there are no specialisations.
	702
	703	\ingroup blas3
	704
	705	\internal
	706
	707	template parameters:
	708	\param M type of the result matrix \c M
	709	\param E1 type of a matrix expression \c A
	710	\param E2 type of a matrix expression \c X
	711	*/
	712	template<class M, class E1, class E2>
	713	BOOST_UBLAS_INLINE
	714	M &
	715	axpy_prod (const matrix_expression<E1> &e1,
	716	const matrix_expression<E2> &e2,
	717	M &m, bool init = true) {
	718	typedef typename M::value_type value_type;
	719	typedef typename M::storage_category storage_category;
	720	typedef typename M::orientation_category orientation_category;
	721
	722	if (init)
	723	m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
	724	return axpy_prod (e1, e2, m, full (), storage_category (), orientation_category ());
	725	}
	726	template<class M, class E1, class E2>
	727	BOOST_UBLAS_INLINE
	728	M
	729	axpy_prod (const matrix_expression<E1> &e1,
	730	const matrix_expression<E2> &e2) {
	731	typedef M matrix_type;
	732
	733	matrix_type m (e1 ().size1 (), e2 ().size2 ());
	734	return axpy_prod (e1, e2, m, full (), true);
	735	}
	736
	737
	738	template<class M, class E1, class E2>
	739	BOOST_UBLAS_INLINE
	740	M &
	741	opb_prod (const matrix_expression<E1> &e1,
	742	const matrix_expression<E2> &e2,
	743	M &m,
	744	dense_proxy_tag, row_major_tag) {
	745	typedef M matrix_type;
	746	typedef const E1 expression1_type;
	747	typedef const E2 expression2_type;
	748	typedef typename M::size_type size_type;
	749	typedef typename M::value_type value_type;
	750
	751	#if BOOST_UBLAS_TYPE_CHECK
	752	matrix<value_type, row_major> cm (m);
	753	typedef typename type_traits<value_type>::real_type real_type;
	754	real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
	755	indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
	756	#endif
	757	size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
	758	for (size_type k = 0; k < size; ++ k) {
	759	vector<value_type> ce1 (column (e1 (), k));
	760	vector<value_type> re2 (row (e2 (), k));
	761	m.plus_assign (outer_prod (ce1, re2));
	762	}
	763	#if BOOST_UBLAS_TYPE_CHECK
	764	BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
	765	#endif
	766	return m;
	767	}
	768
	769	template<class M, class E1, class E2>
	770	BOOST_UBLAS_INLINE
	771	M &
	772	opb_prod (const matrix_expression<E1> &e1,
	773	const matrix_expression<E2> &e2,
	774	M &m,
	775	dense_proxy_tag, column_major_tag) {
	776	typedef M matrix_type;
	777	typedef const E1 expression1_type;
	778	typedef const E2 expression2_type;
	779	typedef typename M::size_type size_type;
	780	typedef typename M::value_type value_type;
	781
	782	#if BOOST_UBLAS_TYPE_CHECK
	783	matrix<value_type, column_major> cm (m);
	784	typedef typename type_traits<value_type>::real_type real_type;
	785	real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
	786	indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
	787	#endif
	788	size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
	789	for (size_type k = 0; k < size; ++ k) {
	790	vector<value_type> ce1 (column (e1 (), k));
	791	vector<value_type> re2 (row (e2 (), k));
	792	m.plus_assign (outer_prod (ce1, re2));
	793	}
	794	#if BOOST_UBLAS_TYPE_CHECK
	795	BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
	796	#endif
	797	return m;
	798	}
	799
	800	// Dispatcher
	801
	802	/** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
	803	optimized fashion.
	804
	805	\param e1 the matrix expression \c A
	806	\param e2 the matrix expression \c X
	807	\param m the result matrix \c M
	808	\param init a boolean parameter
	809
	810	<tt>opb_prod(A, X, M, init)</tt> implements the well known
	811	axpy-product. Setting \a init to \c true is equivalent to call
	812	<tt>M.clear()</tt> before <tt>opb_prod</tt>. Currently \a init
	813	defaults to \c true, but this may change in the future.
	814
	815	This function may give a speedup if \c A has less columns than
	816	rows, because the product is computed as a sum of outer
	817	products.
	818
	819	\ingroup blas3
	820
	821	\internal
	822
	823	template parameters:
	824	\param M type of the result matrix \c M
	825	\param E1 type of a matrix expression \c A
	826	\param E2 type of a matrix expression \c X
	827	*/
	828	template<class M, class E1, class E2>
	829	BOOST_UBLAS_INLINE
	830	M &
	831	opb_prod (const matrix_expression<E1> &e1,
	832	const matrix_expression<E2> &e2,
	833	M &m, bool init = true) {
	834	typedef typename M::value_type value_type;
	835	typedef typename M::storage_category storage_category;
	836	typedef typename M::orientation_category orientation_category;
	837
	838	if (init)
	839	m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
	840	return opb_prod (e1, e2, m, storage_category (), orientation_category ());
	841	}
	842	template<class M, class E1, class E2>
	843	BOOST_UBLAS_INLINE
	844	M
	845	opb_prod (const matrix_expression<E1> &e1,
	846	const matrix_expression<E2> &e2) {
	847	typedef M matrix_type;
	848
	849	matrix_type m (e1 ().size1 (), e2 ().size2 ());
	850	return opb_prod (e1, e2, m, true);
	851	}
	852
	853	}}}
	854
	855	#endif

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

Context Navigation

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

Download in other formats: