Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

operation.hpp @ 857

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

Line
1	//
2	// Copyright (c) 2000-2002
3	// Joerg Walter, Mathias Koch
4	//
5	// Permission to use, copy, modify, distribute and sell this software
6	// and its documentation for any purpose is hereby granted without fee,
7	// provided that the above copyright notice appear in all copies and
8	// that both that copyright notice and this permission notice appear
9	// in supporting documentation. The authors make no representations
10	// about the suitability of this software for any purpose.
11	// It is provided "as is" without express or implied warranty.
12	//
13	// The authors gratefully acknowledge the support of
14	// GeNeSys mbH & Co. KG in producing this work.
15	//
16
17	#ifndef _BOOST_UBLAS_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: