2 // Copyright (c) 2000-2002
3 // Joerg Walter, Mathias Koch
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)
9 // The authors gratefully acknowledge the support of
10 // GeNeSys mbH & Co. KG in producing this work.
13 #ifndef _BOOST_UBLAS_MATRIX_ASSIGN_
14 #define _BOOST_UBLAS_MATRIX_ASSIGN_
16 #include <boost/numeric/ublas/traits.hpp>
17 // Required for make_conformant storage
20 // Iterators based on ideas of Jeremy Siek
22 namespace boost { namespace numeric { namespace ublas {
25 // Weak equality check - useful to compare equality two arbitary matrix expression results.
26 // Since the actual expressions are unknown, we check for and arbitary error bound
27 // on the relative error.
28 // For a linear expression the infinity norm makes sense as we do not know how the elements will be
29 // combined in the expression. False positive results are inevitable for arbirary expressions!
30 template<class E1, class E2, class S>
32 bool equals (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2, S epsilon, S min_norm) {
33 return norm_inf (e1 - e2) < epsilon *
34 std::max<S> (std::max<S> (norm_inf (e1), norm_inf (e2)), min_norm);
37 template<class E1, class E2>
39 bool expression_type_check (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2) {
40 typedef typename type_traits<typename promote_traits<typename E1::value_type,
41 typename E2::value_type>::promote_type>::real_type real_type;
42 return equals (e1, e2, BOOST_UBLAS_TYPE_CHECK_EPSILON, BOOST_UBLAS_TYPE_CHECK_MIN);
46 template<class M, class E, class R>
47 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
48 void make_conformant (M &m, const matrix_expression<E> &e, row_major_tag, R) {
49 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
50 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
51 typedef R conformant_restrict_type;
52 typedef typename M::size_type size_type;
53 typedef typename M::difference_type difference_type;
54 typedef typename M::value_type value_type;
55 // FIXME unbounded_array with push_back maybe better
56 std::vector<std::pair<size_type, size_type> > index;
57 typename M::iterator1 it1 (m.begin1 ());
58 typename M::iterator1 it1_end (m.end1 ());
59 typename E::const_iterator1 it1e (e ().begin1 ());
60 typename E::const_iterator1 it1e_end (e ().end1 ());
61 while (it1 != it1_end && it1e != it1e_end) {
62 difference_type compare = it1.index1 () - it1e.index1 ();
64 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
65 typename M::iterator2 it2 (it1.begin ());
66 typename M::iterator2 it2_end (it1.end ());
67 typename E::const_iterator2 it2e (it1e.begin ());
68 typename E::const_iterator2 it2e_end (it1e.end ());
70 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
71 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
72 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
73 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
75 if (it2 != it2_end && it2e != it2e_end) {
76 size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
78 difference_type compare2 = it2_index - it2e_index;
81 if (it2 != it2_end && it2e != it2e_end) {
82 it2_index = it2.index2 ();
83 it2e_index = it2e.index2 ();
86 } else if (compare2 < 0) {
87 increment (it2, it2_end, - compare2);
89 it2_index = it2.index2 ();
92 } else if (compare2 > 0) {
93 if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
94 if (static_cast<value_type>(*it2e) != value_type/*zero*/())
95 index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
98 it2e_index = it2e.index2 ();
104 while (it2e != it2e_end) {
105 if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
106 if (static_cast<value_type>(*it2e) != value_type/*zero*/())
107 index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
111 } else if (compare < 0) {
112 increment (it1, it1_end, - compare);
113 } else if (compare > 0) {
114 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
115 typename E::const_iterator2 it2e (it1e.begin ());
116 typename E::const_iterator2 it2e_end (it1e.end ());
118 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
119 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
121 while (it2e != it2e_end) {
122 if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
123 if (static_cast<value_type>(*it2e) != value_type/*zero*/())
124 index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
130 while (it1e != it1e_end) {
131 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
132 typename E::const_iterator2 it2e (it1e.begin ());
133 typename E::const_iterator2 it2e_end (it1e.end ());
135 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
136 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
138 while (it2e != it2e_end) {
139 if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
140 if (static_cast<value_type>(*it2e) != value_type/*zero*/())
141 index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
146 // ISSUE proxies require insert_element
147 for (size_type k = 0; k < index.size (); ++ k)
148 m (index [k].first, index [k].second) = value_type/*zero*/();
150 template<class M, class E, class R>
151 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
152 void make_conformant (M &m, const matrix_expression<E> &e, column_major_tag, R) {
153 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
154 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
155 typedef R conformant_restrict_type;
156 typedef typename M::size_type size_type;
157 typedef typename M::difference_type difference_type;
158 typedef typename M::value_type value_type;
159 std::vector<std::pair<size_type, size_type> > index;
160 typename M::iterator2 it2 (m.begin2 ());
161 typename M::iterator2 it2_end (m.end2 ());
162 typename E::const_iterator2 it2e (e ().begin2 ());
163 typename E::const_iterator2 it2e_end (e ().end2 ());
164 while (it2 != it2_end && it2e != it2e_end) {
165 difference_type compare = it2.index2 () - it2e.index2 ();
167 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
168 typename M::iterator1 it1 (it2.begin ());
169 typename M::iterator1 it1_end (it2.end ());
170 typename E::const_iterator1 it1e (it2e.begin ());
171 typename E::const_iterator1 it1e_end (it2e.end ());
173 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
174 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
175 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
176 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
178 if (it1 != it1_end && it1e != it1e_end) {
179 size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
181 difference_type compare2 = it1_index - it1e_index;
184 if (it1 != it1_end && it1e != it1e_end) {
185 it1_index = it1.index1 ();
186 it1e_index = it1e.index1 ();
189 } else if (compare2 < 0) {
190 increment (it1, it1_end, - compare2);
192 it1_index = it1.index1 ();
195 } else if (compare2 > 0) {
196 if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
197 if (static_cast<value_type>(*it1e) != value_type/*zero*/())
198 index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
200 if (it1e != it1e_end)
201 it1e_index = it1e.index1 ();
207 while (it1e != it1e_end) {
208 if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
209 if (static_cast<value_type>(*it1e) != value_type/*zero*/())
210 index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
214 } else if (compare < 0) {
215 increment (it2, it2_end, - compare);
216 } else if (compare > 0) {
217 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
218 typename E::const_iterator1 it1e (it2e.begin ());
219 typename E::const_iterator1 it1e_end (it2e.end ());
221 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
222 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
224 while (it1e != it1e_end) {
225 if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
226 if (static_cast<value_type>(*it1e) != value_type/*zero*/())
227 index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
233 while (it2e != it2e_end) {
234 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
235 typename E::const_iterator1 it1e (it2e.begin ());
236 typename E::const_iterator1 it1e_end (it2e.end ());
238 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
239 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
241 while (it1e != it1e_end) {
242 if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
243 if (static_cast<value_type>(*it1e) != value_type/*zero*/())
244 index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
249 // ISSUE proxies require insert_element
250 for (size_type k = 0; k < index.size (); ++ k)
251 m (index [k].first, index [k].second) = value_type/*zero*/();
257 // Explicitly iterating row major
258 template<template <class T1, class T2> class F, class M, class T>
259 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
260 void iterating_matrix_assign_scalar (M &m, const T &t, row_major_tag) {
261 typedef F<typename M::iterator2::reference, T> functor_type;
262 typedef typename M::difference_type difference_type;
263 difference_type size1 (m.size1 ());
264 difference_type size2 (m.size2 ());
265 typename M::iterator1 it1 (m.begin1 ());
266 BOOST_UBLAS_CHECK (size2 == 0 || m.end1 () - it1 == size1, bad_size ());
267 while (-- size1 >= 0) {
268 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
269 typename M::iterator2 it2 (it1.begin ());
271 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
273 BOOST_UBLAS_CHECK (it1.end () - it2 == size2, bad_size ());
274 difference_type temp_size2 (size2);
275 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
276 while (-- temp_size2 >= 0)
277 functor_type::apply (*it2, t), ++ it2;
279 DD (temp_size2, 4, r, (functor_type::apply (*it2, t), ++ it2));
284 // Explicitly iterating column major
285 template<template <class T1, class T2> class F, class M, class T>
286 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
287 void iterating_matrix_assign_scalar (M &m, const T &t, column_major_tag) {
288 typedef F<typename M::iterator1::reference, T> functor_type;
289 typedef typename M::difference_type difference_type;
290 difference_type size2 (m.size2 ());
291 difference_type size1 (m.size1 ());
292 typename M::iterator2 it2 (m.begin2 ());
293 BOOST_UBLAS_CHECK (size1 == 0 || m.end2 () - it2 == size2, bad_size ());
294 while (-- size2 >= 0) {
295 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
296 typename M::iterator1 it1 (it2.begin ());
298 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
300 BOOST_UBLAS_CHECK (it2.end () - it1 == size1, bad_size ());
301 difference_type temp_size1 (size1);
302 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
303 while (-- temp_size1 >= 0)
304 functor_type::apply (*it1, t), ++ it1;
306 DD (temp_size1, 4, r, (functor_type::apply (*it1, t), ++ it1));
311 // Explicitly indexing row major
312 template<template <class T1, class T2> class F, class M, class T>
313 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
314 void indexing_matrix_assign_scalar (M &m, const T &t, row_major_tag) {
315 typedef F<typename M::reference, T> functor_type;
316 typedef typename M::size_type size_type;
317 size_type size1 (m.size1 ());
318 size_type size2 (m.size2 ());
319 for (size_type i = 0; i < size1; ++ i) {
320 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
321 for (size_type j = 0; j < size2; ++ j)
322 functor_type::apply (m (i, j), t);
325 DD (size2, 4, r, (functor_type::apply (m (i, j), t), ++ j));
329 // Explicitly indexing column major
330 template<template <class T1, class T2> class F, class M, class T>
331 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
332 void indexing_matrix_assign_scalar (M &m, const T &t, column_major_tag) {
333 typedef F<typename M::reference, T> functor_type;
334 typedef typename M::size_type size_type;
335 size_type size2 (m.size2 ());
336 size_type size1 (m.size1 ());
337 for (size_type j = 0; j < size2; ++ j) {
338 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
339 for (size_type i = 0; i < size1; ++ i)
340 functor_type::apply (m (i, j), t);
343 DD (size1, 4, r, (functor_type::apply (m (i, j), t), ++ i));
348 // Dense (proxy) case
349 template<template <class T1, class T2> class F, class M, class T, class C>
350 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
351 void matrix_assign_scalar (M &m, const T &t, dense_proxy_tag, C) {
352 typedef C orientation_category;
353 #ifdef BOOST_UBLAS_USE_INDEXING
354 indexing_matrix_assign_scalar<F> (m, t, orientation_category ());
355 #elif BOOST_UBLAS_USE_ITERATING
356 iterating_matrix_assign_scalar<F> (m, t, orientation_category ());
358 typedef typename M::size_type size_type;
359 size_type size1 (m.size1 ());
360 size_type size2 (m.size2 ());
361 if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
362 size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
363 iterating_matrix_assign_scalar<F> (m, t, orientation_category ());
365 indexing_matrix_assign_scalar<F> (m, t, orientation_category ());
368 // Packed (proxy) row major case
369 template<template <class T1, class T2> class F, class M, class T>
370 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
371 void matrix_assign_scalar (M &m, const T &t, packed_proxy_tag, row_major_tag) {
372 typedef F<typename M::iterator2::reference, T> functor_type;
373 typedef typename M::difference_type difference_type;
374 typename M::iterator1 it1 (m.begin1 ());
375 difference_type size1 (m.end1 () - it1);
376 while (-- size1 >= 0) {
377 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
378 typename M::iterator2 it2 (it1.begin ());
379 difference_type size2 (it1.end () - it2);
381 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
382 difference_type size2 (end (it1, iterator1_tag ()) - it2);
384 while (-- size2 >= 0)
385 functor_type::apply (*it2, t), ++ it2;
389 // Packed (proxy) column major case
390 template<template <class T1, class T2> class F, class M, class T>
391 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
392 void matrix_assign_scalar (M &m, const T &t, packed_proxy_tag, column_major_tag) {
393 typedef F<typename M::iterator1::reference, T> functor_type;
394 typedef typename M::difference_type difference_type;
395 typename M::iterator2 it2 (m.begin2 ());
396 difference_type size2 (m.end2 () - it2);
397 while (-- size2 >= 0) {
398 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
399 typename M::iterator1 it1 (it2.begin ());
400 difference_type size1 (it2.end () - it1);
402 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
403 difference_type size1 (end (it2, iterator2_tag ()) - it1);
405 while (-- size1 >= 0)
406 functor_type::apply (*it1, t), ++ it1;
410 // Sparse (proxy) row major case
411 template<template <class T1, class T2> class F, class M, class T>
412 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
413 void matrix_assign_scalar (M &m, const T &t, sparse_proxy_tag, row_major_tag) {
414 typedef F<typename M::iterator2::reference, T> functor_type;
415 typename M::iterator1 it1 (m.begin1 ());
416 typename M::iterator1 it1_end (m.end1 ());
417 while (it1 != it1_end) {
418 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
419 typename M::iterator2 it2 (it1.begin ());
420 typename M::iterator2 it2_end (it1.end ());
422 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
423 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
425 while (it2 != it2_end)
426 functor_type::apply (*it2, t), ++ it2;
430 // Sparse (proxy) column major case
431 template<template <class T1, class T2> class F, class M, class T>
432 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
433 void matrix_assign_scalar (M &m, const T &t, sparse_proxy_tag, column_major_tag) {
434 typedef F<typename M::iterator1::reference, T> functor_type;
435 typename M::iterator2 it2 (m.begin2 ());
436 typename M::iterator2 it2_end (m.end2 ());
437 while (it2 != it2_end) {
438 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
439 typename M::iterator1 it1 (it2.begin ());
440 typename M::iterator1 it1_end (it2.end ());
442 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
443 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
445 while (it1 != it1_end)
446 functor_type::apply (*it1, t), ++ it1;
452 template<template <class T1, class T2> class F, class M, class T>
454 void matrix_assign_scalar (M &m, const T &t) {
455 typedef typename M::storage_category storage_category;
456 typedef typename M::orientation_category orientation_category;
457 matrix_assign_scalar<F> (m, t, storage_category (), orientation_category ());
460 template<class SC, bool COMPUTED, class RI1, class RI2>
461 struct matrix_assign_traits {
462 typedef SC storage_category;
465 template<bool COMPUTED>
466 struct matrix_assign_traits<dense_tag, COMPUTED, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
467 typedef packed_tag storage_category;
470 struct matrix_assign_traits<dense_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
471 typedef sparse_tag storage_category;
474 struct matrix_assign_traits<dense_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
475 typedef sparse_proxy_tag storage_category;
478 template<bool COMPUTED>
479 struct matrix_assign_traits<dense_proxy_tag, COMPUTED, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
480 typedef packed_proxy_tag storage_category;
482 template<bool COMPUTED>
483 struct matrix_assign_traits<dense_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
484 typedef sparse_proxy_tag storage_category;
488 struct matrix_assign_traits<packed_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
489 typedef sparse_tag storage_category;
492 struct matrix_assign_traits<packed_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
493 typedef sparse_proxy_tag storage_category;
496 template<bool COMPUTED>
497 struct matrix_assign_traits<packed_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
498 typedef sparse_proxy_tag storage_category;
502 struct matrix_assign_traits<sparse_tag, true, dense_random_access_iterator_tag, dense_random_access_iterator_tag> {
503 typedef sparse_proxy_tag storage_category;
506 struct matrix_assign_traits<sparse_tag, true, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
507 typedef sparse_proxy_tag storage_category;
510 struct matrix_assign_traits<sparse_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
511 typedef sparse_proxy_tag storage_category;
514 // Explicitly iterating row major
515 template<template <class T1, class T2> class F, class M, class E>
516 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
517 void iterating_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
518 typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
519 typedef typename M::difference_type difference_type;
520 difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
521 difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
522 typename M::iterator1 it1 (m.begin1 ());
523 BOOST_UBLAS_CHECK (size2 == 0 || m.end1 () - it1 == size1, bad_size ());
524 typename E::const_iterator1 it1e (e ().begin1 ());
525 BOOST_UBLAS_CHECK (size2 == 0 || e ().end1 () - it1e == size1, bad_size ());
526 while (-- size1 >= 0) {
527 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
528 typename M::iterator2 it2 (it1.begin ());
529 typename E::const_iterator2 it2e (it1e.begin ());
531 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
532 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
534 BOOST_UBLAS_CHECK (it1.end () - it2 == size2, bad_size ());
535 BOOST_UBLAS_CHECK (it1e.end () - it2e == size2, bad_size ());
536 difference_type temp_size2 (size2);
537 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
538 while (-- temp_size2 >= 0)
539 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
541 DD (temp_size2, 2, r, (functor_type::apply (*it2, *it2e), ++ it2, ++ it2e));
546 // Explicitly iterating column major
547 template<template <class T1, class T2> class F, class M, class E>
548 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
549 void iterating_matrix_assign (M &m, const matrix_expression<E> &e, column_major_tag) {
550 typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
551 typedef typename M::difference_type difference_type;
552 difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
553 difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
554 typename M::iterator2 it2 (m.begin2 ());
555 BOOST_UBLAS_CHECK (size1 == 0 || m.end2 () - it2 == size2, bad_size ());
556 typename E::const_iterator2 it2e (e ().begin2 ());
557 BOOST_UBLAS_CHECK (size1 == 0 || e ().end2 () - it2e == size2, bad_size ());
558 while (-- size2 >= 0) {
559 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
560 typename M::iterator1 it1 (it2.begin ());
561 typename E::const_iterator1 it1e (it2e.begin ());
563 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
564 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
566 BOOST_UBLAS_CHECK (it2.end () - it1 == size1, bad_size ());
567 BOOST_UBLAS_CHECK (it2e.end () - it1e == size1, bad_size ());
568 difference_type temp_size1 (size1);
569 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
570 while (-- temp_size1 >= 0)
571 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
573 DD (temp_size1, 2, r, (functor_type::apply (*it1, *it1e), ++ it1, ++ it1e));
578 // Explicitly indexing row major
579 template<template <class T1, class T2> class F, class M, class E>
580 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
581 void indexing_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
582 typedef F<typename M::reference, typename E::value_type> functor_type;
583 typedef typename M::size_type size_type;
584 size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
585 size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
586 for (size_type i = 0; i < size1; ++ i) {
587 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
588 for (size_type j = 0; j < size2; ++ j)
589 functor_type::apply (m (i, j), e () (i, j));
592 DD (size2, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ j));
596 // Explicitly indexing column major
597 template<template <class T1, class T2> class F, class M, class E>
598 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
599 void indexing_matrix_assign (M &m, const matrix_expression<E> &e, column_major_tag) {
600 typedef F<typename M::reference, typename E::value_type> functor_type;
601 typedef typename M::size_type size_type;
602 size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
603 size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
604 for (size_type j = 0; j < size2; ++ j) {
605 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
606 for (size_type i = 0; i < size1; ++ i)
607 functor_type::apply (m (i, j), e () (i, j));
610 DD (size1, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ i));
615 // Dense (proxy) case
616 template<template <class T1, class T2> class F, class R, class M, class E, class C>
617 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
618 void matrix_assign (M &m, const matrix_expression<E> &e, dense_proxy_tag, C) {
619 // R unnecessary, make_conformant not required
620 typedef C orientation_category;
621 #ifdef BOOST_UBLAS_USE_INDEXING
622 indexing_matrix_assign<F> (m, e, orientation_category ());
623 #elif BOOST_UBLAS_USE_ITERATING
624 iterating_matrix_assign<F> (m, e, orientation_category ());
626 typedef typename M::difference_type difference_type;
627 size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
628 size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
629 if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
630 size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
631 iterating_matrix_assign<F> (m, e, orientation_category ());
633 indexing_matrix_assign<F> (m, e, orientation_category ());
636 // Packed (proxy) row major case
637 template<template <class T1, class T2> class F, class R, class M, class E>
638 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
639 void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
640 typedef typename matrix_traits<E>::value_type expr_value_type;
641 typedef F<typename M::iterator2::reference, expr_value_type> functor_type;
642 // R unnecessary, make_conformant not required
643 typedef typename M::difference_type difference_type;
645 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
646 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
648 #if BOOST_UBLAS_TYPE_CHECK
649 typedef typename M::value_type value_type;
650 matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
651 indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ());
652 indexing_matrix_assign<F> (cm, e, row_major_tag ());
654 typename M::iterator1 it1 (m.begin1 ());
655 typename M::iterator1 it1_end (m.end1 ());
656 typename E::const_iterator1 it1e (e ().begin1 ());
657 typename E::const_iterator1 it1e_end (e ().end1 ());
658 difference_type it1_size (it1_end - it1);
659 difference_type it1e_size (it1e_end - it1e);
660 difference_type diff1 (0);
661 if (it1_size > 0 && it1e_size > 0)
662 diff1 = it1.index1 () - it1e.index1 ();
664 difference_type size1 = (std::min) (diff1, it1e_size);
670 size1 = (std::min) (- diff1, it1_size);
673 if (!functor_type::computed) {
674 while (-- size1 >= 0) { // zeroing
675 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
676 typename M::iterator2 it2 (it1.begin ());
677 typename M::iterator2 it2_end (it1.end ());
679 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
680 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
682 difference_type size2 (it2_end - it2);
683 while (-- size2 >= 0)
684 functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
693 difference_type size1 ((std::min) (it1_size, it1e_size));
696 while (-- size1 >= 0) {
697 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
698 typename M::iterator2 it2 (it1.begin ());
699 typename M::iterator2 it2_end (it1.end ());
700 typename E::const_iterator2 it2e (it1e.begin ());
701 typename E::const_iterator2 it2e_end (it1e.end ());
703 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
704 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
705 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
706 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
708 difference_type it2_size (it2_end - it2);
709 difference_type it2e_size (it2e_end - it2e);
710 difference_type diff2 (0);
711 if (it2_size > 0 && it2e_size > 0) {
712 diff2 = it2.index2 () - it2e.index2 ();
713 difference_type size2 = (std::min) (diff2, it2e_size);
719 size2 = (std::min) (- diff2, it2_size);
722 if (!functor_type::computed) {
723 while (-- size2 >= 0) // zeroing
724 functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
731 difference_type size2 ((std::min) (it2_size, it2e_size));
734 while (-- size2 >= 0)
735 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
737 if (!functor_type::computed) {
738 while (-- size2 >= 0) // zeroing
739 functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
746 if (!functor_type::computed) {
747 while (-- size1 >= 0) { // zeroing
748 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
749 typename M::iterator2 it2 (it1.begin ());
750 typename M::iterator2 it2_end (it1.end ());
752 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
753 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
755 difference_type size2 (it2_end - it2);
756 while (-- size2 >= 0)
757 functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
763 #if BOOST_UBLAS_TYPE_CHECK
764 if (! disable_type_check<bool>::value)
765 BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
768 // Packed (proxy) column major case
769 template<template <class T1, class T2> class F, class R, class M, class E>
770 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
771 void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, column_major_tag) {
772 typedef typename matrix_traits<E>::value_type expr_value_type;
773 typedef F<typename M::iterator1::reference, expr_value_type> functor_type;
774 // R unnecessary, make_conformant not required
775 typedef typename M::difference_type difference_type;
777 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
778 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
780 #if BOOST_UBLAS_TYPE_CHECK
781 typedef typename M::value_type value_type;
782 matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
783 indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ());
784 indexing_matrix_assign<F> (cm, e, column_major_tag ());
786 typename M::iterator2 it2 (m.begin2 ());
787 typename M::iterator2 it2_end (m.end2 ());
788 typename E::const_iterator2 it2e (e ().begin2 ());
789 typename E::const_iterator2 it2e_end (e ().end2 ());
790 difference_type it2_size (it2_end - it2);
791 difference_type it2e_size (it2e_end - it2e);
792 difference_type diff2 (0);
793 if (it2_size > 0 && it2e_size > 0)
794 diff2 = it2.index2 () - it2e.index2 ();
796 difference_type size2 = (std::min) (diff2, it2e_size);
802 size2 = (std::min) (- diff2, it2_size);
805 if (!functor_type::computed) {
806 while (-- size2 >= 0) { // zeroing
807 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
808 typename M::iterator1 it1 (it2.begin ());
809 typename M::iterator1 it1_end (it2.end ());
811 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
812 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
814 difference_type size1 (it1_end - it1);
815 while (-- size1 >= 0)
816 functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
825 difference_type size2 ((std::min) (it2_size, it2e_size));
828 while (-- size2 >= 0) {
829 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
830 typename M::iterator1 it1 (it2.begin ());
831 typename M::iterator1 it1_end (it2.end ());
832 typename E::const_iterator1 it1e (it2e.begin ());
833 typename E::const_iterator1 it1e_end (it2e.end ());
835 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
836 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
837 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
838 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
840 difference_type it1_size (it1_end - it1);
841 difference_type it1e_size (it1e_end - it1e);
842 difference_type diff1 (0);
843 if (it1_size > 0 && it1e_size > 0) {
844 diff1 = it1.index1 () - it1e.index1 ();
845 difference_type size1 = (std::min) (diff1, it1e_size);
851 size1 = (std::min) (- diff1, it1_size);
854 if (!functor_type::computed) {
855 while (-- size1 >= 0) // zeroing
856 functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
863 difference_type size1 ((std::min) (it1_size, it1e_size));
866 while (-- size1 >= 0)
867 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
869 if (!functor_type::computed) {
870 while (-- size1 >= 0) // zeroing
871 functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
878 if (!functor_type::computed) {
879 while (-- size2 >= 0) { // zeroing
880 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
881 typename M::iterator1 it1 (it2.begin ());
882 typename M::iterator1 it1_end (it2.end ());
884 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
885 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
887 difference_type size1 (it1_end - it1);
888 while (-- size1 >= 0)
889 functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
895 #if BOOST_UBLAS_TYPE_CHECK
896 if (! disable_type_check<bool>::value)
897 BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
900 // Sparse row major case
901 template<template <class T1, class T2> class F, class R, class M, class E>
902 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
903 void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, row_major_tag) {
904 typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
905 // R unnecessary, make_conformant not required
906 BOOST_STATIC_ASSERT ((!functor_type::computed));
907 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
908 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
909 typedef typename M::value_type value_type;
910 // Sparse type has no numeric constraints to check
913 typename E::const_iterator1 it1e (e ().begin1 ());
914 typename E::const_iterator1 it1e_end (e ().end1 ());
915 while (it1e != it1e_end) {
916 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
917 typename E::const_iterator2 it2e (it1e.begin ());
918 typename E::const_iterator2 it2e_end (it1e.end ());
920 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
921 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
923 while (it2e != it2e_end) {
924 value_type t (*it2e);
925 if (t != value_type/*zero*/())
926 m.insert_element (it2e.index1 (), it2e.index2 (), t);
932 // Sparse column major case
933 template<template <class T1, class T2> class F, class R, class M, class E>
934 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
935 void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, column_major_tag) {
936 typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
937 // R unnecessary, make_conformant not required
938 BOOST_STATIC_ASSERT ((!functor_type::computed));
939 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
940 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
941 typedef typename M::value_type value_type;
942 // Sparse type has no numeric constraints to check
945 typename E::const_iterator2 it2e (e ().begin2 ());
946 typename E::const_iterator2 it2e_end (e ().end2 ());
947 while (it2e != it2e_end) {
948 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
949 typename E::const_iterator1 it1e (it2e.begin ());
950 typename E::const_iterator1 it1e_end (it2e.end ());
952 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
953 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
955 while (it1e != it1e_end) {
956 value_type t (*it1e);
957 if (t != value_type/*zero*/())
958 m.insert_element (it1e.index1 (), it1e.index2 (), t);
964 // Sparse proxy or functional row major case
965 template<template <class T1, class T2> class F, class R, class M, class E>
966 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
967 void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) {
968 typedef typename matrix_traits<E>::value_type expr_value_type;
969 typedef F<typename M::iterator2::reference, expr_value_type> functor_type;
970 typedef R conformant_restrict_type;
971 typedef typename M::size_type size_type;
972 typedef typename M::difference_type difference_type;
974 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
975 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
977 #if BOOST_UBLAS_TYPE_CHECK
978 typedef typename M::value_type value_type;
979 matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
980 indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ());
981 indexing_matrix_assign<F> (cm, e, row_major_tag ());
983 detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
985 typename M::iterator1 it1 (m.begin1 ());
986 typename M::iterator1 it1_end (m.end1 ());
987 typename E::const_iterator1 it1e (e ().begin1 ());
988 typename E::const_iterator1 it1e_end (e ().end1 ());
989 while (it1 != it1_end && it1e != it1e_end) {
990 difference_type compare = it1.index1 () - it1e.index1 ();
992 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
993 typename M::iterator2 it2 (it1.begin ());
994 typename M::iterator2 it2_end (it1.end ());
995 typename E::const_iterator2 it2e (it1e.begin ());
996 typename E::const_iterator2 it2e_end (it1e.end ());
998 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
999 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1000 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
1001 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
1003 if (it2 != it2_end && it2e != it2e_end) {
1004 size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
1006 difference_type compare2 = it2_index - it2e_index;
1007 if (compare2 == 0) {
1008 functor_type::apply (*it2, *it2e);
1010 if (it2 != it2_end && it2e != it2e_end) {
1011 it2_index = it2.index2 ();
1012 it2e_index = it2e.index2 ();
1015 } else if (compare2 < 0) {
1016 if (!functor_type::computed) {
1017 functor_type::apply (*it2, expr_value_type/*zero*/());
1020 increment (it2, it2_end, - compare2);
1022 it2_index = it2.index2 ();
1025 } else if (compare2 > 0) {
1026 increment (it2e, it2e_end, compare2);
1027 if (it2e != it2e_end)
1028 it2e_index = it2e.index2 ();
1034 if (!functor_type::computed) {
1035 while (it2 != it2_end) { // zeroing
1036 functor_type::apply (*it2, expr_value_type/*zero*/());
1043 } else if (compare < 0) {
1044 if (!functor_type::computed) {
1045 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1046 typename M::iterator2 it2 (it1.begin ());
1047 typename M::iterator2 it2_end (it1.end ());
1049 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1050 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1052 while (it2 != it2_end) { // zeroing
1053 functor_type::apply (*it2, expr_value_type/*zero*/());
1058 increment (it1, it1_end, - compare);
1060 } else if (compare > 0) {
1061 increment (it1e, it1e_end, compare);
1064 if (!functor_type::computed) {
1065 while (it1 != it1_end) {
1066 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1067 typename M::iterator2 it2 (it1.begin ());
1068 typename M::iterator2 it2_end (it1.end ());
1070 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1071 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1073 while (it2 != it2_end) { // zeroing
1074 functor_type::apply (*it2, expr_value_type/*zero*/());
1082 #if BOOST_UBLAS_TYPE_CHECK
1083 if (! disable_type_check<bool>::value)
1084 BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
1087 // Sparse proxy or functional column major case
1088 template<template <class T1, class T2> class F, class R, class M, class E>
1089 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1090 void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) {
1091 typedef typename matrix_traits<E>::value_type expr_value_type;
1092 typedef F<typename M::iterator1::reference, expr_value_type> functor_type;
1093 typedef R conformant_restrict_type;
1094 typedef typename M::size_type size_type;
1095 typedef typename M::difference_type difference_type;
1097 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
1098 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
1100 #if BOOST_UBLAS_TYPE_CHECK
1101 typedef typename M::value_type value_type;
1102 matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
1103 indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ());
1104 indexing_matrix_assign<F> (cm, e, column_major_tag ());
1106 detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
1108 typename M::iterator2 it2 (m.begin2 ());
1109 typename M::iterator2 it2_end (m.end2 ());
1110 typename E::const_iterator2 it2e (e ().begin2 ());
1111 typename E::const_iterator2 it2e_end (e ().end2 ());
1112 while (it2 != it2_end && it2e != it2e_end) {
1113 difference_type compare = it2.index2 () - it2e.index2 ();
1115 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1116 typename M::iterator1 it1 (it2.begin ());
1117 typename M::iterator1 it1_end (it2.end ());
1118 typename E::const_iterator1 it1e (it2e.begin ());
1119 typename E::const_iterator1 it1e_end (it2e.end ());
1121 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1122 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1123 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
1124 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
1126 if (it1 != it1_end && it1e != it1e_end) {
1127 size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
1129 difference_type compare2 = it1_index - it1e_index;
1130 if (compare2 == 0) {
1131 functor_type::apply (*it1, *it1e);
1133 if (it1 != it1_end && it1e != it1e_end) {
1134 it1_index = it1.index1 ();
1135 it1e_index = it1e.index1 ();
1138 } else if (compare2 < 0) {
1139 if (!functor_type::computed) {
1140 functor_type::apply (*it1, expr_value_type/*zero*/()); // zeroing
1143 increment (it1, it1_end, - compare2);
1145 it1_index = it1.index1 ();
1148 } else if (compare2 > 0) {
1149 increment (it1e, it1e_end, compare2);
1150 if (it1e != it1e_end)
1151 it1e_index = it1e.index1 ();
1157 if (!functor_type::computed) {
1158 while (it1 != it1_end) { // zeroing
1159 functor_type::apply (*it1, expr_value_type/*zero*/());
1166 } else if (compare < 0) {
1167 if (!functor_type::computed) {
1168 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1169 typename M::iterator1 it1 (it2.begin ());
1170 typename M::iterator1 it1_end (it2.end ());
1172 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1173 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1175 while (it1 != it1_end) { // zeroing
1176 functor_type::apply (*it1, expr_value_type/*zero*/());
1181 increment (it2, it2_end, - compare);
1183 } else if (compare > 0) {
1184 increment (it2e, it2e_end, compare);
1187 if (!functor_type::computed) {
1188 while (it2 != it2_end) {
1189 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1190 typename M::iterator1 it1 (it2.begin ());
1191 typename M::iterator1 it1_end (it2.end ());
1193 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1194 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1196 while (it1 != it1_end) { // zeroing
1197 functor_type::apply (*it1, expr_value_type/*zero*/());
1205 #if BOOST_UBLAS_TYPE_CHECK
1206 if (! disable_type_check<bool>::value)
1207 BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
1212 template<template <class T1, class T2> class F, class M, class E>
1214 void matrix_assign (M &m, const matrix_expression<E> &e) {
1215 typedef typename matrix_assign_traits<typename M::storage_category,
1216 F<typename M::reference, typename E::value_type>::computed,
1217 typename E::const_iterator1::iterator_category,
1218 typename E::const_iterator2::iterator_category>::storage_category storage_category;
1219 // give preference to matrix M's orientation if known
1220 typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
1221 typename E::orientation_category ,
1222 typename M::orientation_category >::type orientation_category;
1223 typedef basic_full<typename M::size_type> unrestricted;
1224 matrix_assign<F, unrestricted> (m, e, storage_category (), orientation_category ());
1226 template<template <class T1, class T2> class F, class R, class M, class E>
1228 void matrix_assign (M &m, const matrix_expression<E> &e) {
1229 typedef R conformant_restrict_type;
1230 typedef typename matrix_assign_traits<typename M::storage_category,
1231 F<typename M::reference, typename E::value_type>::computed,
1232 typename E::const_iterator1::iterator_category,
1233 typename E::const_iterator2::iterator_category>::storage_category storage_category;
1234 // give preference to matrix M's orientation if known
1235 typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
1236 typename E::orientation_category ,
1237 typename M::orientation_category >::type orientation_category;
1238 matrix_assign<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ());
1241 template<class SC, class RI1, class RI2>
1242 struct matrix_swap_traits {
1243 typedef SC storage_category;
1247 struct matrix_swap_traits<dense_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
1248 typedef sparse_proxy_tag storage_category;
1252 struct matrix_swap_traits<packed_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
1253 typedef sparse_proxy_tag storage_category;
1256 // Dense (proxy) row major case
1257 template<template <class T1, class T2> class F, class R, class M, class E>
1258 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1259 void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, row_major_tag) {
1260 typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
1261 // R unnecessary, make_conformant not required
1262 //typedef typename M::size_type size_type; // gcc is complaining that this is not used, although this is not right
1263 typedef typename M::difference_type difference_type;
1264 typename M::iterator1 it1 (m.begin1 ());
1265 typename E::iterator1 it1e (e ().begin1 ());
1266 difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (e ().end1 () - it1e)));
1267 while (-- size1 >= 0) {
1268 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1269 typename M::iterator2 it2 (it1.begin ());
1270 typename E::iterator2 it2e (it1e.begin ());
1271 difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (it1e.end () - it2e)));
1273 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1274 typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1275 difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (end (it1e, iterator1_tag ()) - it2e)));
1277 while (-- size2 >= 0)
1278 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
1282 // Dense (proxy) column major case
1283 template<template <class T1, class T2> class F, class R, class M, class E>
1284 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1285 void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, column_major_tag) {
1286 typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
1287 // R unnecessary, make_conformant not required
1288 // typedef typename M::size_type size_type; // gcc is complaining that this is not used, although this is not right
1289 typedef typename M::difference_type difference_type;
1290 typename M::iterator2 it2 (m.begin2 ());
1291 typename E::iterator2 it2e (e ().begin2 ());
1292 difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (e ().end2 () - it2e)));
1293 while (-- size2 >= 0) {
1294 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1295 typename M::iterator1 it1 (it2.begin ());
1296 typename E::iterator1 it1e (it2e.begin ());
1297 difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (it2e.end () - it1e)));
1299 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1300 typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1301 difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (end (it2e, iterator2_tag ()) - it1e)));
1303 while (-- size1 >= 0)
1304 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
1308 // Packed (proxy) row major case
1309 template<template <class T1, class T2> class F, class R, class M, class E>
1310 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1311 void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
1312 typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
1313 // R unnecessary, make_conformant not required
1314 typedef typename M::difference_type difference_type;
1315 typename M::iterator1 it1 (m.begin1 ());
1316 typename E::iterator1 it1e (e ().begin1 ());
1317 difference_type size1 (BOOST_UBLAS_SAME (m.end1 () - it1, e ().end1 () - it1e));
1318 while (-- size1 >= 0) {
1319 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1320 typename M::iterator2 it2 (it1.begin ());
1321 typename E::iterator2 it2e (it1e.begin ());
1322 difference_type size2 (BOOST_UBLAS_SAME (it1.end () - it2, it1e.end () - it2e));
1324 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1325 typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1326 difference_type size2 (BOOST_UBLAS_SAME (end (it1, iterator1_tag ()) - it2, end (it1e, iterator1_tag ()) - it2e));
1328 while (-- size2 >= 0)
1329 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
1333 // Packed (proxy) column major case
1334 template<template <class T1, class T2> class F, class R, class M, class E>
1335 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1336 void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, column_major_tag) {
1337 typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
1338 // R unnecessary, make_conformant not required
1339 typedef typename M::difference_type difference_type;
1340 typename M::iterator2 it2 (m.begin2 ());
1341 typename E::iterator2 it2e (e ().begin2 ());
1342 difference_type size2 (BOOST_UBLAS_SAME (m.end2 () - it2, e ().end2 () - it2e));
1343 while (-- size2 >= 0) {
1344 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1345 typename M::iterator1 it1 (it2.begin ());
1346 typename E::iterator1 it1e (it2e.begin ());
1347 difference_type size1 (BOOST_UBLAS_SAME (it2.end () - it1, it2e.end () - it1e));
1349 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1350 typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1351 difference_type size1 (BOOST_UBLAS_SAME (end (it2, iterator2_tag ()) - it1, end (it2e, iterator2_tag ()) - it1e));
1353 while (-- size1 >= 0)
1354 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
1358 // Sparse (proxy) row major case
1359 template<template <class T1, class T2> class F, class R, class M, class E>
1360 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1361 void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) {
1362 typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
1363 typedef R conformant_restrict_type;
1364 typedef typename M::size_type size_type;
1365 typedef typename M::difference_type difference_type;
1366 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
1367 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
1369 detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
1370 // FIXME should be a seperate restriction for E
1371 detail::make_conformant (e (), m, row_major_tag (), conformant_restrict_type ());
1373 typename M::iterator1 it1 (m.begin1 ());
1374 typename M::iterator1 it1_end (m.end1 ());
1375 typename E::iterator1 it1e (e ().begin1 ());
1376 typename E::iterator1 it1e_end (e ().end1 ());
1377 while (it1 != it1_end && it1e != it1e_end) {
1378 difference_type compare = it1.index1 () - it1e.index1 ();
1380 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1381 typename M::iterator2 it2 (it1.begin ());
1382 typename M::iterator2 it2_end (it1.end ());
1383 typename E::iterator2 it2e (it1e.begin ());
1384 typename E::iterator2 it2e_end (it1e.end ());
1386 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1387 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1388 typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1389 typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
1391 if (it2 != it2_end && it2e != it2e_end) {
1392 size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
1394 difference_type compare2 = it2_index - it2e_index;
1395 if (compare2 == 0) {
1396 functor_type::apply (*it2, *it2e);
1398 if (it2 != it2_end && it2e != it2e_end) {
1399 it2_index = it2.index2 ();
1400 it2e_index = it2e.index2 ();
1403 } else if (compare2 < 0) {
1404 increment (it2, it2_end, - compare2);
1406 it2_index = it2.index2 ();
1409 } else if (compare2 > 0) {
1410 increment (it2e, it2e_end, compare2);
1411 if (it2e != it2e_end)
1412 it2e_index = it2e.index2 ();
1418 #if BOOST_UBLAS_TYPE_CHECK
1419 increment (it2e, it2e_end);
1420 increment (it2, it2_end);
1423 } else if (compare < 0) {
1424 #if BOOST_UBLAS_TYPE_CHECK
1425 while (it1.index1 () < it1e.index1 ()) {
1426 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1427 typename M::iterator2 it2 (it1.begin ());
1428 typename M::iterator2 it2_end (it1.end ());
1430 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1431 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1433 increment (it2, it2_end);
1437 increment (it1, it1_end, - compare);
1439 } else if (compare > 0) {
1440 #if BOOST_UBLAS_TYPE_CHECK
1441 while (it1e.index1 () < it1.index1 ()) {
1442 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1443 typename E::iterator2 it2e (it1e.begin ());
1444 typename E::iterator2 it2e_end (it1e.end ());
1446 typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1447 typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
1449 increment (it2e, it2e_end);
1453 increment (it1e, it1e_end, compare);
1457 #if BOOST_UBLAS_TYPE_CHECK
1458 while (it1e != it1e_end) {
1459 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1460 typename E::iterator2 it2e (it1e.begin ());
1461 typename E::iterator2 it2e_end (it1e.end ());
1463 typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1464 typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
1466 increment (it2e, it2e_end);
1469 while (it1 != it1_end) {
1470 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1471 typename M::iterator2 it2 (it1.begin ());
1472 typename M::iterator2 it2_end (it1.end ());
1474 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1475 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1477 increment (it2, it2_end);
1482 // Sparse (proxy) column major case
1483 template<template <class T1, class T2> class F, class R, class M, class E>
1484 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1485 void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) {
1486 typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
1487 typedef R conformant_restrict_type;
1488 typedef typename M::size_type size_type;
1489 typedef typename M::difference_type difference_type;
1491 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
1492 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
1494 detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
1495 // FIXME should be a seperate restriction for E
1496 detail::make_conformant (e (), m, column_major_tag (), conformant_restrict_type ());
1498 typename M::iterator2 it2 (m.begin2 ());
1499 typename M::iterator2 it2_end (m.end2 ());
1500 typename E::iterator2 it2e (e ().begin2 ());
1501 typename E::iterator2 it2e_end (e ().end2 ());
1502 while (it2 != it2_end && it2e != it2e_end) {
1503 difference_type compare = it2.index2 () - it2e.index2 ();
1505 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1506 typename M::iterator1 it1 (it2.begin ());
1507 typename M::iterator1 it1_end (it2.end ());
1508 typename E::iterator1 it1e (it2e.begin ());
1509 typename E::iterator1 it1e_end (it2e.end ());
1511 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1512 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1513 typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1514 typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
1516 if (it1 != it1_end && it1e != it1e_end) {
1517 size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
1519 difference_type compare2 = it1_index - it1e_index;
1520 if (compare2 == 0) {
1521 functor_type::apply (*it1, *it1e);
1523 if (it1 != it1_end && it1e != it1e_end) {
1524 it1_index = it1.index1 ();
1525 it1e_index = it1e.index1 ();
1528 } else if (compare2 < 0) {
1529 increment (it1, it1_end, - compare2);
1531 it1_index = it1.index1 ();
1534 } else if (compare2 > 0) {
1535 increment (it1e, it1e_end, compare2);
1536 if (it1e != it1e_end)
1537 it1e_index = it1e.index1 ();
1543 #if BOOST_UBLAS_TYPE_CHECK
1544 increment (it1e, it1e_end);
1545 increment (it1, it1_end);
1548 } else if (compare < 0) {
1549 #if BOOST_UBLAS_TYPE_CHECK
1550 while (it2.index2 () < it2e.index2 ()) {
1551 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1552 typename M::iterator1 it1 (it2.begin ());
1553 typename M::iterator1 it1_end (it2.end ());
1555 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1556 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1558 increment (it1, it1_end);
1562 increment (it2, it2_end, - compare);
1564 } else if (compare > 0) {
1565 #if BOOST_UBLAS_TYPE_CHECK
1566 while (it2e.index2 () < it2.index2 ()) {
1567 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1568 typename E::iterator1 it1e (it2e.begin ());
1569 typename E::iterator1 it1e_end (it2e.end ());
1571 typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1572 typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
1574 increment (it1e, it1e_end);
1578 increment (it2e, it2e_end, compare);
1582 #if BOOST_UBLAS_TYPE_CHECK
1583 while (it2e != it2e_end) {
1584 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1585 typename E::iterator1 it1e (it2e.begin ());
1586 typename E::iterator1 it1e_end (it2e.end ());
1588 typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1589 typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
1591 increment (it1e, it1e_end);
1594 while (it2 != it2_end) {
1595 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1596 typename M::iterator1 it1 (it2.begin ());
1597 typename M::iterator1 it1_end (it2.end ());
1599 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1600 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1602 increment (it1, it1_end);
1609 template<template <class T1, class T2> class F, class M, class E>
1611 void matrix_swap (M &m, matrix_expression<E> &e) {
1612 typedef typename matrix_swap_traits<typename M::storage_category,
1613 typename E::const_iterator1::iterator_category,
1614 typename E::const_iterator2::iterator_category>::storage_category storage_category;
1615 // give preference to matrix M's orientation if known
1616 typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
1617 typename E::orientation_category ,
1618 typename M::orientation_category >::type orientation_category;
1619 typedef basic_full<typename M::size_type> unrestricted;
1620 matrix_swap<F, unrestricted> (m, e, storage_category (), orientation_category ());
1622 template<template <class T1, class T2> class F, class R, class M, class E>
1624 void matrix_swap (M &m, matrix_expression<E> &e) {
1625 typedef R conformant_restrict_type;
1626 typedef typename matrix_swap_traits<typename M::storage_category,
1627 typename E::const_iterator1::iterator_category,
1628 typename E::const_iterator2::iterator_category>::storage_category storage_category;
1629 // give preference to matrix M's orientation if known
1630 typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
1631 typename E::orientation_category ,
1632 typename M::orientation_category >::type orientation_category;
1633 matrix_swap<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ());