]>
Commit | Line | Data |
---|---|---|
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 | ||
26 | namespace 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 |