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

source code of include/boost/numeric/ublas/operation.hpp