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_MATRIX_ASSIGN_ 14 #define _BOOST_UBLAS_MATRIX_ASSIGN_ 15 16 #include <boost/numeric/ublas/traits.hpp> 17 // Required for make_conformant storage 18 #include <vector> 19 20 // Iterators based on ideas of Jeremy Siek 21 22 namespace boost { namespace numeric { namespace ublas { 23 namespace detail { 24 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> 31 BOOST_UBLAS_INLINE equals(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,S epsilon,S min_norm)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); 35 } 36 37 template<class E1, class E2> 38 BOOST_UBLAS_INLINE expression_type_check(const matrix_expression<E1> & e1,const matrix_expression<E2> & 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); 43 } 44 45 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. make_conformant(M & m,const matrix_expression<E> & e,row_major_tag,R)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 (); 63 if (compare == 0) { 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 ()); 69 #else 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 ())); 74 #endif 75 if (it2 != it2_end && it2e != it2e_end) { 76 size_type it2_index = it2.index2 (), it2e_index = it2e.index2 (); 77 for (;;) { 78 difference_type compare2 = it2_index - it2e_index; 79 if (compare2 == 0) { 80 ++ it2, ++ it2e; 81 if (it2 != it2_end && it2e != it2e_end) { 82 it2_index = it2.index2 (); 83 it2e_index = it2e.index2 (); 84 } else 85 break; 86 } else if (compare2 < 0) { 87 increment (it2, it2_end, - compare2); 88 if (it2 != it2_end) 89 it2_index = it2.index2 (); 90 else 91 break; 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 ())); 96 ++ it2e; 97 if (it2e != it2e_end) 98 it2e_index = it2e.index2 (); 99 else 100 break; 101 } 102 } 103 } 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 ())); 108 ++ it2e; 109 } 110 ++ it1, ++ it1e; 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 ()); 117 #else 118 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); 119 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ())); 120 #endif 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 ())); 125 ++ it2e; 126 } 127 ++ it1e; 128 } 129 } 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 ()); 134 #else 135 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); 136 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ())); 137 #endif 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 ())); 142 ++ it2e; 143 } 144 ++ it1e; 145 } 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*/(); 149 } 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. make_conformant(M & m,const matrix_expression<E> & e,column_major_tag,R)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 (); 166 if (compare == 0) { 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 ()); 172 #else 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 ())); 177 #endif 178 if (it1 != it1_end && it1e != it1e_end) { 179 size_type it1_index = it1.index1 (), it1e_index = it1e.index1 (); 180 for (;;) { 181 difference_type compare2 = it1_index - it1e_index; 182 if (compare2 == 0) { 183 ++ it1, ++ it1e; 184 if (it1 != it1_end && it1e != it1e_end) { 185 it1_index = it1.index1 (); 186 it1e_index = it1e.index1 (); 187 } else 188 break; 189 } else if (compare2 < 0) { 190 increment (it1, it1_end, - compare2); 191 if (it1 != it1_end) 192 it1_index = it1.index1 (); 193 else 194 break; 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 ())); 199 ++ it1e; 200 if (it1e != it1e_end) 201 it1e_index = it1e.index1 (); 202 else 203 break; 204 } 205 } 206 } 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 ())); 211 ++ it1e; 212 } 213 ++ it2, ++ it2e; 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 ()); 220 #else 221 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); 222 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ())); 223 #endif 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 ())); 228 ++ it1e; 229 } 230 ++ it2e; 231 } 232 } 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 ()); 237 #else 238 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); 239 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ())); 240 #endif 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 ())); 245 ++ it1e; 246 } 247 ++ it2e; 248 } 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*/(); 252 } 253 254 }//namespace detail 255 256 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. iterating_matrix_assign_scalar(M & m,const T & t,row_major_tag)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 ()); 270 #else 271 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 272 #endif 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; 278 #else 279 DD (temp_size2, 4, r, (functor_type::apply (*it2, t), ++ it2)); 280 #endif 281 ++ it1; 282 } 283 } 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. iterating_matrix_assign_scalar(M & m,const T & t,column_major_tag)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 ()); 297 #else 298 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 299 #endif 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; 305 #else 306 DD (temp_size1, 4, r, (functor_type::apply (*it1, t), ++ it1)); 307 #endif 308 ++ it2; 309 } 310 } 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. indexing_matrix_assign_scalar(M & m,const T & t,row_major_tag)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); 323 #else 324 size_type j (0); 325 DD (size2, 4, r, (functor_type::apply (m (i, j), t), ++ j)); 326 #endif 327 } 328 } 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. indexing_matrix_assign_scalar(M & m,const T & t,column_major_tag)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); 341 #else 342 size_type i (0); 343 DD (size1, 4, r, (functor_type::apply (m (i, j), t), ++ i)); 344 #endif 345 } 346 } 347 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. matrix_assign_scalar(M & m,const T & t,dense_proxy_tag,C)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 ()); 357 #else 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 ()); 364 else 365 indexing_matrix_assign_scalar<F> (m, t, orientation_category ()); 366 #endif 367 } 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. matrix_assign_scalar(M & m,const T & t,packed_proxy_tag,row_major_tag)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); 380 #else 381 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 382 difference_type size2 (end (it1, iterator1_tag ()) - it2); 383 #endif 384 while (-- size2 >= 0) 385 functor_type::apply (*it2, t), ++ it2; 386 ++ it1; 387 } 388 } 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. matrix_assign_scalar(M & m,const T & t,packed_proxy_tag,column_major_tag)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); 401 #else 402 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 403 difference_type size1 (end (it2, iterator2_tag ()) - it1); 404 #endif 405 while (-- size1 >= 0) 406 functor_type::apply (*it1, t), ++ it1; 407 ++ it2; 408 } 409 } 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. matrix_assign_scalar(M & m,const T & t,sparse_proxy_tag,row_major_tag)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 ()); 421 #else 422 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 423 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 424 #endif 425 while (it2 != it2_end) 426 functor_type::apply (*it2, t), ++ it2; 427 ++ it1; 428 } 429 } 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. matrix_assign_scalar(M & m,const T & t,sparse_proxy_tag,column_major_tag)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 ()); 441 #else 442 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 443 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 444 #endif 445 while (it1 != it1_end) 446 functor_type::apply (*it1, t), ++ it1; 447 ++ it2; 448 } 449 } 450 451 // Dispatcher 452 template<template <class T1, class T2> class F, class M, class T> 453 BOOST_UBLAS_INLINE matrix_assign_scalar(M & m,const T & 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 ()); 458 } 459 460 template<class SC, bool COMPUTED, class RI1, class RI2> 461 struct matrix_assign_traits { 462 typedef SC storage_category; 463 }; 464 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; 468 }; 469 template<> 470 struct matrix_assign_traits<dense_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> { 471 typedef sparse_tag storage_category; 472 }; 473 template<> 474 struct matrix_assign_traits<dense_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> { 475 typedef sparse_proxy_tag storage_category; 476 }; 477 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; 481 }; 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; 485 }; 486 487 template<> 488 struct matrix_assign_traits<packed_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> { 489 typedef sparse_tag storage_category; 490 }; 491 template<> 492 struct matrix_assign_traits<packed_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> { 493 typedef sparse_proxy_tag storage_category; 494 }; 495 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; 499 }; 500 501 template<> 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; 504 }; 505 template<> 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; 508 }; 509 template<> 510 struct matrix_assign_traits<sparse_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> { 511 typedef sparse_proxy_tag storage_category; 512 }; 513 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. iterating_matrix_assign(M & m,const matrix_expression<E> & e,row_major_tag)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 ()); 530 #else 531 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 532 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); 533 #endif 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; 540 #else 541 DD (temp_size2, 2, r, (functor_type::apply (*it2, *it2e), ++ it2, ++ it2e)); 542 #endif 543 ++ it1, ++ it1e; 544 } 545 } 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. iterating_matrix_assign(M & m,const matrix_expression<E> & e,column_major_tag)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 ()); 562 #else 563 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 564 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); 565 #endif 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; 572 #else 573 DD (temp_size1, 2, r, (functor_type::apply (*it1, *it1e), ++ it1, ++ it1e)); 574 #endif 575 ++ it2, ++ it2e; 576 } 577 } 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. indexing_matrix_assign(M & m,const matrix_expression<E> & e,row_major_tag)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)); 590 #else 591 size_type j (0); 592 DD (size2, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ j)); 593 #endif 594 } 595 } 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. indexing_matrix_assign(M & m,const matrix_expression<E> & e,column_major_tag)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)); 608 #else 609 size_type i (0); 610 DD (size1, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ i)); 611 #endif 612 } 613 } 614 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. matrix_assign(M & m,const matrix_expression<E> & e,dense_proxy_tag,C)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 ()); 625 #else 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 ()); 632 else 633 indexing_matrix_assign<F> (m, e, orientation_category ()); 634 #endif 635 } 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. matrix_assign(M & m,const matrix_expression<E> & e,packed_proxy_tag,row_major_tag)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; 644 645 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); 646 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); 647 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 ()); 653 #endif 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 (); 663 if (diff1 != 0) { 664 difference_type size1 = (std::min) (diff1, it1e_size); 665 if (size1 > 0) { 666 it1e += size1; 667 it1e_size -= size1; 668 diff1 -= size1; 669 } 670 size1 = (std::min) (- diff1, it1_size); 671 if (size1 > 0) { 672 it1_size -= size1; 673 //Disabled warning C4127 because the conditional expression is constant 674 #ifdef _MSC_VER 675 #pragma warning(push) 676 #pragma warning(disable: 4127) 677 #endif 678 if (!functor_type::computed) { 679 #ifdef _MSC_VER 680 #pragma warning(pop) 681 #endif 682 while (-- size1 >= 0) { // zeroing 683 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 684 typename M::iterator2 it2 (it1.begin ()); 685 typename M::iterator2 it2_end (it1.end ()); 686 #else 687 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 688 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 689 #endif 690 difference_type size2 (it2_end - it2); 691 while (-- size2 >= 0) 692 functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2; 693 ++ it1; 694 } 695 } else { 696 it1 += size1; 697 } 698 diff1 += size1; 699 } 700 } 701 difference_type size1 ((std::min) (it1_size, it1e_size)); 702 it1_size -= size1; 703 it1e_size -= size1; 704 while (-- size1 >= 0) { 705 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 706 typename M::iterator2 it2 (it1.begin ()); 707 typename M::iterator2 it2_end (it1.end ()); 708 typename E::const_iterator2 it2e (it1e.begin ()); 709 typename E::const_iterator2 it2e_end (it1e.end ()); 710 #else 711 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 712 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 713 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); 714 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ())); 715 #endif 716 difference_type it2_size (it2_end - it2); 717 difference_type it2e_size (it2e_end - it2e); 718 difference_type diff2 (0); 719 if (it2_size > 0 && it2e_size > 0) { 720 diff2 = it2.index2 () - it2e.index2 (); 721 difference_type size2 = (std::min) (diff2, it2e_size); 722 if (size2 > 0) { 723 it2e += size2; 724 it2e_size -= size2; 725 diff2 -= size2; 726 } 727 size2 = (std::min) (- diff2, it2_size); 728 if (size2 > 0) { 729 it2_size -= size2; 730 //Disabled warning C4127 because the conditional expression is constant 731 #ifdef _MSC_VER 732 #pragma warning(push) 733 #pragma warning(disable: 4127) 734 #endif 735 if (!functor_type::computed) { 736 #ifdef _MSC_VER 737 #pragma warning(pop) 738 #endif 739 while (-- size2 >= 0) // zeroing 740 functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2; 741 } else { 742 it2 += size2; 743 } 744 diff2 += size2; 745 } 746 } 747 difference_type size2 ((std::min) (it2_size, it2e_size)); 748 it2_size -= size2; 749 it2e_size -= size2; 750 while (-- size2 >= 0) 751 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e; 752 size2 = it2_size; 753 //Disabled warning C4127 because the conditional expression is constant 754 #ifdef _MSC_VER 755 #pragma warning(push) 756 #pragma warning(disable: 4127) 757 #endif 758 if (!functor_type::computed) { 759 #ifdef _MSC_VER 760 #pragma warning(pop) 761 #endif 762 while (-- size2 >= 0) // zeroing 763 functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2; 764 } else { 765 it2 += size2; 766 } 767 ++ it1, ++ it1e; 768 } 769 size1 = it1_size; 770 //Disabled warning C4127 because the conditional expression is constant 771 #ifdef _MSC_VER 772 #pragma warning(push) 773 #pragma warning(disable: 4127) 774 #endif 775 if (!functor_type::computed) { 776 #ifdef _MSC_VER 777 #pragma warning(pop) 778 #endif 779 while (-- size1 >= 0) { // zeroing 780 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 781 typename M::iterator2 it2 (it1.begin ()); 782 typename M::iterator2 it2_end (it1.end ()); 783 #else 784 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 785 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 786 #endif 787 difference_type size2 (it2_end - it2); 788 while (-- size2 >= 0) 789 functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2; 790 ++ it1; 791 } 792 } else { 793 it1 += size1; 794 } 795 #if BOOST_UBLAS_TYPE_CHECK 796 if (! disable_type_check<bool>::value) 797 BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ()); 798 #endif 799 } 800 // Packed (proxy) column major case 801 template<template <class T1, class T2> class F, class R, class M, class E> 802 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_assign(M & m,const matrix_expression<E> & e,packed_proxy_tag,column_major_tag)803 void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, column_major_tag) { 804 typedef typename matrix_traits<E>::value_type expr_value_type; 805 typedef F<typename M::iterator1::reference, expr_value_type> functor_type; 806 // R unnecessary, make_conformant not required 807 typedef typename M::difference_type difference_type; 808 809 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); 810 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); 811 812 #if BOOST_UBLAS_TYPE_CHECK 813 typedef typename M::value_type value_type; 814 matrix<value_type, column_major> cm (m.size1 (), m.size2 ()); 815 indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ()); 816 indexing_matrix_assign<F> (cm, e, column_major_tag ()); 817 #endif 818 typename M::iterator2 it2 (m.begin2 ()); 819 typename M::iterator2 it2_end (m.end2 ()); 820 typename E::const_iterator2 it2e (e ().begin2 ()); 821 typename E::const_iterator2 it2e_end (e ().end2 ()); 822 difference_type it2_size (it2_end - it2); 823 difference_type it2e_size (it2e_end - it2e); 824 difference_type diff2 (0); 825 if (it2_size > 0 && it2e_size > 0) 826 diff2 = it2.index2 () - it2e.index2 (); 827 if (diff2 != 0) { 828 difference_type size2 = (std::min) (diff2, it2e_size); 829 if (size2 > 0) { 830 it2e += size2; 831 it2e_size -= size2; 832 diff2 -= size2; 833 } 834 size2 = (std::min) (- diff2, it2_size); 835 if (size2 > 0) { 836 it2_size -= size2; 837 //Disabled warning C4127 because the conditional expression is constant 838 #ifdef _MSC_VER 839 #pragma warning(push) 840 #pragma warning(disable: 4127) 841 #endif 842 if (!functor_type::computed) { 843 #ifdef _MSC_VER 844 #pragma warning(pop) 845 #endif 846 while (-- size2 >= 0) { // zeroing 847 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 848 typename M::iterator1 it1 (it2.begin ()); 849 typename M::iterator1 it1_end (it2.end ()); 850 #else 851 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 852 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 853 #endif 854 difference_type size1 (it1_end - it1); 855 while (-- size1 >= 0) 856 functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1; 857 ++ it2; 858 } 859 } else { 860 it2 += size2; 861 } 862 diff2 += size2; 863 } 864 } 865 difference_type size2 ((std::min) (it2_size, it2e_size)); 866 it2_size -= size2; 867 it2e_size -= size2; 868 while (-- size2 >= 0) { 869 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 870 typename M::iterator1 it1 (it2.begin ()); 871 typename M::iterator1 it1_end (it2.end ()); 872 typename E::const_iterator1 it1e (it2e.begin ()); 873 typename E::const_iterator1 it1e_end (it2e.end ()); 874 #else 875 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 876 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 877 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); 878 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ())); 879 #endif 880 difference_type it1_size (it1_end - it1); 881 difference_type it1e_size (it1e_end - it1e); 882 difference_type diff1 (0); 883 if (it1_size > 0 && it1e_size > 0) { 884 diff1 = it1.index1 () - it1e.index1 (); 885 difference_type size1 = (std::min) (diff1, it1e_size); 886 if (size1 > 0) { 887 it1e += size1; 888 it1e_size -= size1; 889 diff1 -= size1; 890 } 891 size1 = (std::min) (- diff1, it1_size); 892 if (size1 > 0) { 893 it1_size -= size1; 894 //Disabled warning C4127 because the conditional expression is constant 895 #ifdef _MSC_VER 896 #pragma warning(push) 897 #pragma warning(disable: 4127) 898 #endif 899 if (!functor_type::computed) { 900 #ifdef _MSC_VER 901 #pragma warning(pop) 902 #endif 903 while (-- size1 >= 0) // zeroing 904 functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1; 905 } else { 906 it1 += size1; 907 } 908 diff1 += size1; 909 } 910 } 911 difference_type size1 ((std::min) (it1_size, it1e_size)); 912 it1_size -= size1; 913 it1e_size -= size1; 914 while (-- size1 >= 0) 915 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e; 916 size1 = it1_size; 917 //Disabled warning C4127 because the conditional expression is constant 918 #ifdef _MSC_VER 919 #pragma warning(push) 920 #pragma warning(disable: 4127) 921 #endif 922 if (!functor_type::computed) { 923 924 #ifdef _MSC_VER 925 #pragma warning(pop) 926 #endif 927 while (-- size1 >= 0) // zeroing 928 functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1; 929 } else { 930 it1 += size1; 931 } 932 ++ it2, ++ it2e; 933 } 934 size2 = it2_size; 935 //Disabled warning C4127 because the conditional expression is constant 936 #ifdef _MSC_VER 937 #pragma warning(push) 938 #pragma warning(disable: 4127) 939 #endif 940 if (!functor_type::computed) { 941 #ifdef _MSC_VER 942 #pragma warning(pop) 943 #endif 944 while (-- size2 >= 0) { // zeroing 945 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 946 typename M::iterator1 it1 (it2.begin ()); 947 typename M::iterator1 it1_end (it2.end ()); 948 #else 949 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 950 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 951 #endif 952 difference_type size1 (it1_end - it1); 953 while (-- size1 >= 0) 954 functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1; 955 ++ it2; 956 } 957 } else { 958 it2 += size2; 959 } 960 #if BOOST_UBLAS_TYPE_CHECK 961 if (! disable_type_check<bool>::value) 962 BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ()); 963 #endif 964 } 965 // Sparse row major case 966 template<template <class T1, class T2> class F, class R, class M, class E> 967 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_assign(M & m,const matrix_expression<E> & e,sparse_tag,row_major_tag)968 void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, row_major_tag) { 969 typedef F<typename M::iterator2::reference, typename E::value_type> functor_type; 970 // R unnecessary, make_conformant not required 971 972 //Disabled warning C4127 because the conditional expression is constant 973 #ifdef _MSC_VER 974 #pragma warning(push) 975 #pragma warning(disable: 4127) 976 #endif 977 BOOST_STATIC_ASSERT ((!functor_type::computed)); 978 #ifdef _MSC_VER 979 #pragma warning(pop) 980 #endif 981 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); 982 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); 983 typedef typename M::value_type value_type; 984 // Sparse type has no numeric constraints to check 985 986 m.clear (); 987 typename E::const_iterator1 it1e (e ().begin1 ()); 988 typename E::const_iterator1 it1e_end (e ().end1 ()); 989 while (it1e != it1e_end) { 990 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 991 typename E::const_iterator2 it2e (it1e.begin ()); 992 typename E::const_iterator2 it2e_end (it1e.end ()); 993 #else 994 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); 995 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ())); 996 #endif 997 while (it2e != it2e_end) { 998 value_type t (*it2e); 999 if (t != value_type/*zero*/()) 1000 m.insert_element (it2e.index1 (), it2e.index2 (), t); 1001 ++ it2e; 1002 } 1003 ++ it1e; 1004 } 1005 } 1006 // Sparse column major case 1007 template<template <class T1, class T2> class F, class R, class M, class E> 1008 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_assign(M & m,const matrix_expression<E> & e,sparse_tag,column_major_tag)1009 void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, column_major_tag) { 1010 typedef F<typename M::iterator1::reference, typename E::value_type> functor_type; 1011 // R unnecessary, make_conformant not required 1012 1013 //Disabled warning C4127 because the conditional expression is constant 1014 #ifdef _MSC_VER 1015 #pragma warning(push) 1016 #pragma warning(disable: 4127) 1017 #endif 1018 BOOST_STATIC_ASSERT ((!functor_type::computed)); 1019 #ifdef _MSC_VER 1020 #pragma warning(pop) 1021 #endif 1022 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); 1023 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); 1024 typedef typename M::value_type value_type; 1025 // Sparse type has no numeric constraints to check 1026 1027 m.clear (); 1028 typename E::const_iterator2 it2e (e ().begin2 ()); 1029 typename E::const_iterator2 it2e_end (e ().end2 ()); 1030 while (it2e != it2e_end) { 1031 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1032 typename E::const_iterator1 it1e (it2e.begin ()); 1033 typename E::const_iterator1 it1e_end (it2e.end ()); 1034 #else 1035 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); 1036 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ())); 1037 #endif 1038 while (it1e != it1e_end) { 1039 value_type t (*it1e); 1040 if (t != value_type/*zero*/()) 1041 m.insert_element (it1e.index1 (), it1e.index2 (), t); 1042 ++ it1e; 1043 } 1044 ++ it2e; 1045 } 1046 } 1047 // Sparse proxy or functional row major case 1048 template<template <class T1, class T2> class F, class R, class M, class E> 1049 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_assign(M & m,const matrix_expression<E> & e,sparse_proxy_tag,row_major_tag)1050 void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) { 1051 typedef typename matrix_traits<E>::value_type expr_value_type; 1052 typedef F<typename M::iterator2::reference, expr_value_type> functor_type; 1053 typedef R conformant_restrict_type; 1054 typedef typename M::size_type size_type; 1055 typedef typename M::difference_type difference_type; 1056 1057 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); 1058 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); 1059 1060 #if BOOST_UBLAS_TYPE_CHECK 1061 typedef typename M::value_type value_type; 1062 matrix<value_type, row_major> cm (m.size1 (), m.size2 ()); 1063 indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ()); 1064 indexing_matrix_assign<F> (cm, e, row_major_tag ()); 1065 #endif 1066 detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ()); 1067 1068 typename M::iterator1 it1 (m.begin1 ()); 1069 typename M::iterator1 it1_end (m.end1 ()); 1070 typename E::const_iterator1 it1e (e ().begin1 ()); 1071 typename E::const_iterator1 it1e_end (e ().end1 ()); 1072 while (it1 != it1_end && it1e != it1e_end) { 1073 difference_type compare = it1.index1 () - it1e.index1 (); 1074 if (compare == 0) { 1075 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1076 typename M::iterator2 it2 (it1.begin ()); 1077 typename M::iterator2 it2_end (it1.end ()); 1078 typename E::const_iterator2 it2e (it1e.begin ()); 1079 typename E::const_iterator2 it2e_end (it1e.end ()); 1080 #else 1081 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 1082 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 1083 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ())); 1084 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ())); 1085 #endif 1086 if (it2 != it2_end && it2e != it2e_end) { 1087 size_type it2_index = it2.index2 (), it2e_index = it2e.index2 (); 1088 for (;;) { 1089 difference_type compare2 = it2_index - it2e_index; 1090 if (compare2 == 0) { 1091 functor_type::apply (*it2, *it2e); 1092 ++ it2, ++ it2e; 1093 if (it2 != it2_end && it2e != it2e_end) { 1094 it2_index = it2.index2 (); 1095 it2e_index = it2e.index2 (); 1096 } else 1097 break; 1098 } else if (compare2 < 0) { 1099 //Disabled warning C4127 because the conditional expression is constant 1100 #ifdef _MSC_VER 1101 #pragma warning(push) 1102 #pragma warning(disable: 4127) 1103 #endif 1104 if (!functor_type::computed) { 1105 #ifdef _MSC_VER 1106 #pragma warning(pop) 1107 #endif 1108 functor_type::apply (*it2, expr_value_type/*zero*/()); 1109 ++ it2; 1110 } else 1111 increment (it2, it2_end, - compare2); 1112 if (it2 != it2_end) 1113 it2_index = it2.index2 (); 1114 else 1115 break; 1116 } else if (compare2 > 0) { 1117 increment (it2e, it2e_end, compare2); 1118 if (it2e != it2e_end) 1119 it2e_index = it2e.index2 (); 1120 else 1121 break; 1122 } 1123 } 1124 } 1125 //Disabled warning C4127 because the conditional expression is constant 1126 #ifdef _MSC_VER 1127 #pragma warning(push) 1128 #pragma warning(disable: 4127) 1129 #endif 1130 if (!functor_type::computed) { 1131 #ifdef _MSC_VER 1132 #pragma warning(pop) 1133 #endif 1134 while (it2 != it2_end) { // zeroing 1135 functor_type::apply (*it2, expr_value_type/*zero*/()); 1136 ++ it2; 1137 } 1138 } else { 1139 it2 = it2_end; 1140 } 1141 ++ it1, ++ it1e; 1142 } else if (compare < 0) { 1143 //Disabled warning C4127 because the conditional expression is constant 1144 #ifdef _MSC_VER 1145 #pragma warning(push) 1146 #pragma warning(disable: 4127) 1147 #endif 1148 if (!functor_type::computed) { 1149 #ifdef _MSC_VER 1150 #pragma warning(pop) 1151 #endif 1152 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1153 typename M::iterator2 it2 (it1.begin ()); 1154 typename M::iterator2 it2_end (it1.end ()); 1155 #else 1156 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 1157 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 1158 #endif 1159 while (it2 != it2_end) { // zeroing 1160 functor_type::apply (*it2, expr_value_type/*zero*/()); 1161 ++ it2; 1162 } 1163 ++ it1; 1164 } else { 1165 increment (it1, it1_end, - compare); 1166 } 1167 } else if (compare > 0) { 1168 increment (it1e, it1e_end, compare); 1169 } 1170 } 1171 //Disabled warning C4127 because the conditional expression is constant 1172 #ifdef _MSC_VER 1173 #pragma warning(push) 1174 #pragma warning(disable: 4127) 1175 #endif 1176 if (!functor_type::computed) { 1177 #ifdef _MSC_VER 1178 #pragma warning(pop) 1179 #endif 1180 while (it1 != it1_end) { 1181 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1182 typename M::iterator2 it2 (it1.begin ()); 1183 typename M::iterator2 it2_end (it1.end ()); 1184 #else 1185 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 1186 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 1187 #endif 1188 while (it2 != it2_end) { // zeroing 1189 functor_type::apply (*it2, expr_value_type/*zero*/()); 1190 ++ it2; 1191 } 1192 ++ it1; 1193 } 1194 } else { 1195 it1 = it1_end; 1196 } 1197 #if BOOST_UBLAS_TYPE_CHECK 1198 if (! disable_type_check<bool>::value) 1199 BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ()); 1200 #endif 1201 } 1202 // Sparse proxy or functional column major case 1203 template<template <class T1, class T2> class F, class R, class M, class E> 1204 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_assign(M & m,const matrix_expression<E> & e,sparse_proxy_tag,column_major_tag)1205 void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) { 1206 typedef typename matrix_traits<E>::value_type expr_value_type; 1207 typedef F<typename M::iterator1::reference, expr_value_type> functor_type; 1208 typedef R conformant_restrict_type; 1209 typedef typename M::size_type size_type; 1210 typedef typename M::difference_type difference_type; 1211 1212 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); 1213 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); 1214 1215 #if BOOST_UBLAS_TYPE_CHECK 1216 typedef typename M::value_type value_type; 1217 matrix<value_type, column_major> cm (m.size1 (), m.size2 ()); 1218 indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ()); 1219 indexing_matrix_assign<F> (cm, e, column_major_tag ()); 1220 #endif 1221 detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ()); 1222 1223 typename M::iterator2 it2 (m.begin2 ()); 1224 typename M::iterator2 it2_end (m.end2 ()); 1225 typename E::const_iterator2 it2e (e ().begin2 ()); 1226 typename E::const_iterator2 it2e_end (e ().end2 ()); 1227 while (it2 != it2_end && it2e != it2e_end) { 1228 difference_type compare = it2.index2 () - it2e.index2 (); 1229 if (compare == 0) { 1230 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1231 typename M::iterator1 it1 (it2.begin ()); 1232 typename M::iterator1 it1_end (it2.end ()); 1233 typename E::const_iterator1 it1e (it2e.begin ()); 1234 typename E::const_iterator1 it1e_end (it2e.end ()); 1235 #else 1236 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 1237 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 1238 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ())); 1239 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ())); 1240 #endif 1241 if (it1 != it1_end && it1e != it1e_end) { 1242 size_type it1_index = it1.index1 (), it1e_index = it1e.index1 (); 1243 for (;;) { 1244 difference_type compare2 = it1_index - it1e_index; 1245 if (compare2 == 0) { 1246 functor_type::apply (*it1, *it1e); 1247 ++ it1, ++ it1e; 1248 if (it1 != it1_end && it1e != it1e_end) { 1249 it1_index = it1.index1 (); 1250 it1e_index = it1e.index1 (); 1251 } else 1252 break; 1253 } else if (compare2 < 0) { 1254 //Disabled warning C4127 because the conditional expression is constant 1255 #ifdef _MSC_VER 1256 #pragma warning(push) 1257 #pragma warning(disable: 4127) 1258 #endif 1259 if (!functor_type::computed) { 1260 #ifdef _MSC_VER 1261 #pragma warning(pop) 1262 #endif 1263 functor_type::apply (*it1, expr_value_type/*zero*/()); // zeroing 1264 ++ it1; 1265 } else 1266 increment (it1, it1_end, - compare2); 1267 if (it1 != it1_end) 1268 it1_index = it1.index1 (); 1269 else 1270 break; 1271 } else if (compare2 > 0) { 1272 increment (it1e, it1e_end, compare2); 1273 if (it1e != it1e_end) 1274 it1e_index = it1e.index1 (); 1275 else 1276 break; 1277 } 1278 } 1279 } 1280 //Disabled warning C4127 because the conditional expression is constant 1281 #ifdef _MSC_VER 1282 #pragma warning(push) 1283 #pragma warning(disable: 4127) 1284 #endif 1285 if (!functor_type::computed) { 1286 #ifdef _MSC_VER 1287 #pragma warning(pop) 1288 #endif 1289 while (it1 != it1_end) { // zeroing 1290 functor_type::apply (*it1, expr_value_type/*zero*/()); 1291 ++ it1; 1292 } 1293 } else { 1294 it1 = it1_end; 1295 } 1296 ++ it2, ++ it2e; 1297 } else if (compare < 0) { 1298 //Disabled warning C4127 because the conditional expression is constant 1299 #ifdef _MSC_VER 1300 #pragma warning(push) 1301 #pragma warning(disable: 4127) 1302 #endif 1303 if (!functor_type::computed) { 1304 #ifdef _MSC_VER 1305 #pragma warning(pop) 1306 #endif 1307 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1308 typename M::iterator1 it1 (it2.begin ()); 1309 typename M::iterator1 it1_end (it2.end ()); 1310 #else 1311 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 1312 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 1313 #endif 1314 while (it1 != it1_end) { // zeroing 1315 functor_type::apply (*it1, expr_value_type/*zero*/()); 1316 ++ it1; 1317 } 1318 ++ it2; 1319 } else { 1320 increment (it2, it2_end, - compare); 1321 } 1322 } else if (compare > 0) { 1323 increment (it2e, it2e_end, compare); 1324 } 1325 } 1326 //Disabled warning C4127 because the conditional expression is constant 1327 #ifdef _MSC_VER 1328 #pragma warning(push) 1329 #pragma warning(disable: 4127) 1330 #endif 1331 if (!functor_type::computed) { 1332 #ifdef _MSC_VER 1333 #pragma warning(pop) 1334 #endif 1335 while (it2 != it2_end) { 1336 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1337 typename M::iterator1 it1 (it2.begin ()); 1338 typename M::iterator1 it1_end (it2.end ()); 1339 #else 1340 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 1341 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 1342 #endif 1343 while (it1 != it1_end) { // zeroing 1344 functor_type::apply (*it1, expr_value_type/*zero*/()); 1345 ++ it1; 1346 } 1347 ++ it2; 1348 } 1349 } else { 1350 it2 = it2_end; 1351 } 1352 #if BOOST_UBLAS_TYPE_CHECK 1353 if (! disable_type_check<bool>::value) 1354 BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ()); 1355 #endif 1356 } 1357 1358 // Dispatcher 1359 template<template <class T1, class T2> class F, class M, class E> 1360 BOOST_UBLAS_INLINE matrix_assign(M & m,const matrix_expression<E> & e)1361 void matrix_assign (M &m, const matrix_expression<E> &e) { 1362 typedef typename matrix_assign_traits<typename M::storage_category, 1363 F<typename M::reference, typename E::value_type>::computed, 1364 typename E::const_iterator1::iterator_category, 1365 typename E::const_iterator2::iterator_category>::storage_category storage_category; 1366 // give preference to matrix M's orientation if known 1367 typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>, 1368 typename E::orientation_category , 1369 typename M::orientation_category >::type orientation_category; 1370 typedef basic_full<typename M::size_type> unrestricted; 1371 matrix_assign<F, unrestricted> (m, e, storage_category (), orientation_category ()); 1372 } 1373 template<template <class T1, class T2> class F, class R, class M, class E> 1374 BOOST_UBLAS_INLINE matrix_assign(M & m,const matrix_expression<E> & e)1375 void matrix_assign (M &m, const matrix_expression<E> &e) { 1376 typedef R conformant_restrict_type; 1377 typedef typename matrix_assign_traits<typename M::storage_category, 1378 F<typename M::reference, typename E::value_type>::computed, 1379 typename E::const_iterator1::iterator_category, 1380 typename E::const_iterator2::iterator_category>::storage_category storage_category; 1381 // give preference to matrix M's orientation if known 1382 typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>, 1383 typename E::orientation_category , 1384 typename M::orientation_category >::type orientation_category; 1385 matrix_assign<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ()); 1386 } 1387 1388 template<class SC, class RI1, class RI2> 1389 struct matrix_swap_traits { 1390 typedef SC storage_category; 1391 }; 1392 1393 template<> 1394 struct matrix_swap_traits<dense_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> { 1395 typedef sparse_proxy_tag storage_category; 1396 }; 1397 1398 template<> 1399 struct matrix_swap_traits<packed_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> { 1400 typedef sparse_proxy_tag storage_category; 1401 }; 1402 1403 // Dense (proxy) row major case 1404 template<template <class T1, class T2> class F, class R, class M, class E> 1405 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_swap(M & m,matrix_expression<E> & e,dense_proxy_tag,row_major_tag)1406 void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, row_major_tag) { 1407 typedef F<typename M::iterator2::reference, typename E::reference> functor_type; 1408 // R unnecessary, make_conformant not required 1409 //typedef typename M::size_type size_type; // gcc is complaining that this is not used, although this is not right 1410 typedef typename M::difference_type difference_type; 1411 typename M::iterator1 it1 (m.begin1 ()); 1412 typename E::iterator1 it1e (e ().begin1 ()); 1413 difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (e ().end1 () - it1e))); 1414 while (-- size1 >= 0) { 1415 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1416 typename M::iterator2 it2 (it1.begin ()); 1417 typename E::iterator2 it2e (it1e.begin ()); 1418 difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (it1e.end () - it2e))); 1419 #else 1420 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 1421 typename E::iterator2 it2e (begin (it1e, iterator1_tag ())); 1422 difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (end (it1e, iterator1_tag ()) - it2e))); 1423 #endif 1424 while (-- size2 >= 0) 1425 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e; 1426 ++ it1, ++ it1e; 1427 } 1428 } 1429 // Dense (proxy) column major case 1430 template<template <class T1, class T2> class F, class R, class M, class E> 1431 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_swap(M & m,matrix_expression<E> & e,dense_proxy_tag,column_major_tag)1432 void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, column_major_tag) { 1433 typedef F<typename M::iterator1::reference, typename E::reference> functor_type; 1434 // R unnecessary, make_conformant not required 1435 // typedef typename M::size_type size_type; // gcc is complaining that this is not used, although this is not right 1436 typedef typename M::difference_type difference_type; 1437 typename M::iterator2 it2 (m.begin2 ()); 1438 typename E::iterator2 it2e (e ().begin2 ()); 1439 difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (e ().end2 () - it2e))); 1440 while (-- size2 >= 0) { 1441 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1442 typename M::iterator1 it1 (it2.begin ()); 1443 typename E::iterator1 it1e (it2e.begin ()); 1444 difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (it2e.end () - it1e))); 1445 #else 1446 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 1447 typename E::iterator1 it1e (begin (it2e, iterator2_tag ())); 1448 difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (end (it2e, iterator2_tag ()) - it1e))); 1449 #endif 1450 while (-- size1 >= 0) 1451 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e; 1452 ++ it2, ++ it2e; 1453 } 1454 } 1455 // Packed (proxy) row major case 1456 template<template <class T1, class T2> class F, class R, class M, class E> 1457 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_swap(M & m,matrix_expression<E> & e,packed_proxy_tag,row_major_tag)1458 void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, row_major_tag) { 1459 typedef F<typename M::iterator2::reference, typename E::reference> functor_type; 1460 // R unnecessary, make_conformant not required 1461 typedef typename M::difference_type difference_type; 1462 typename M::iterator1 it1 (m.begin1 ()); 1463 typename E::iterator1 it1e (e ().begin1 ()); 1464 difference_type size1 (BOOST_UBLAS_SAME (m.end1 () - it1, e ().end1 () - it1e)); 1465 while (-- size1 >= 0) { 1466 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1467 typename M::iterator2 it2 (it1.begin ()); 1468 typename E::iterator2 it2e (it1e.begin ()); 1469 difference_type size2 (BOOST_UBLAS_SAME (it1.end () - it2, it1e.end () - it2e)); 1470 #else 1471 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 1472 typename E::iterator2 it2e (begin (it1e, iterator1_tag ())); 1473 difference_type size2 (BOOST_UBLAS_SAME (end (it1, iterator1_tag ()) - it2, end (it1e, iterator1_tag ()) - it2e)); 1474 #endif 1475 while (-- size2 >= 0) 1476 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e; 1477 ++ it1, ++ it1e; 1478 } 1479 } 1480 // Packed (proxy) column major case 1481 template<template <class T1, class T2> class F, class R, class M, class E> 1482 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_swap(M & m,matrix_expression<E> & e,packed_proxy_tag,column_major_tag)1483 void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, column_major_tag) { 1484 typedef F<typename M::iterator1::reference, typename E::reference> functor_type; 1485 // R unnecessary, make_conformant not required 1486 typedef typename M::difference_type difference_type; 1487 typename M::iterator2 it2 (m.begin2 ()); 1488 typename E::iterator2 it2e (e ().begin2 ()); 1489 difference_type size2 (BOOST_UBLAS_SAME (m.end2 () - it2, e ().end2 () - it2e)); 1490 while (-- size2 >= 0) { 1491 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1492 typename M::iterator1 it1 (it2.begin ()); 1493 typename E::iterator1 it1e (it2e.begin ()); 1494 difference_type size1 (BOOST_UBLAS_SAME (it2.end () - it1, it2e.end () - it1e)); 1495 #else 1496 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 1497 typename E::iterator1 it1e (begin (it2e, iterator2_tag ())); 1498 difference_type size1 (BOOST_UBLAS_SAME (end (it2, iterator2_tag ()) - it1, end (it2e, iterator2_tag ()) - it1e)); 1499 #endif 1500 while (-- size1 >= 0) 1501 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e; 1502 ++ it2, ++ it2e; 1503 } 1504 } 1505 // Sparse (proxy) row major case 1506 template<template <class T1, class T2> class F, class R, class M, class E> 1507 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_swap(M & m,matrix_expression<E> & e,sparse_proxy_tag,row_major_tag)1508 void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) { 1509 typedef F<typename M::iterator2::reference, typename E::reference> functor_type; 1510 typedef R conformant_restrict_type; 1511 typedef typename M::size_type size_type; 1512 typedef typename M::difference_type difference_type; 1513 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); 1514 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); 1515 1516 detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ()); 1517 // FIXME should be a seperate restriction for E 1518 detail::make_conformant (e (), m, row_major_tag (), conformant_restrict_type ()); 1519 1520 typename M::iterator1 it1 (m.begin1 ()); 1521 typename M::iterator1 it1_end (m.end1 ()); 1522 typename E::iterator1 it1e (e ().begin1 ()); 1523 typename E::iterator1 it1e_end (e ().end1 ()); 1524 while (it1 != it1_end && it1e != it1e_end) { 1525 difference_type compare = it1.index1 () - it1e.index1 (); 1526 if (compare == 0) { 1527 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1528 typename M::iterator2 it2 (it1.begin ()); 1529 typename M::iterator2 it2_end (it1.end ()); 1530 typename E::iterator2 it2e (it1e.begin ()); 1531 typename E::iterator2 it2e_end (it1e.end ()); 1532 #else 1533 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 1534 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 1535 typename E::iterator2 it2e (begin (it1e, iterator1_tag ())); 1536 typename E::iterator2 it2e_end (end (it1e, iterator1_tag ())); 1537 #endif 1538 if (it2 != it2_end && it2e != it2e_end) { 1539 size_type it2_index = it2.index2 (), it2e_index = it2e.index2 (); 1540 for (;;) { 1541 difference_type compare2 = it2_index - it2e_index; 1542 if (compare2 == 0) { 1543 functor_type::apply (*it2, *it2e); 1544 ++ it2, ++ it2e; 1545 if (it2 != it2_end && it2e != it2e_end) { 1546 it2_index = it2.index2 (); 1547 it2e_index = it2e.index2 (); 1548 } else 1549 break; 1550 } else if (compare2 < 0) { 1551 increment (it2, it2_end, - compare2); 1552 if (it2 != it2_end) 1553 it2_index = it2.index2 (); 1554 else 1555 break; 1556 } else if (compare2 > 0) { 1557 increment (it2e, it2e_end, compare2); 1558 if (it2e != it2e_end) 1559 it2e_index = it2e.index2 (); 1560 else 1561 break; 1562 } 1563 } 1564 } 1565 #if BOOST_UBLAS_TYPE_CHECK 1566 increment (it2e, it2e_end); 1567 increment (it2, it2_end); 1568 #endif 1569 ++ it1, ++ it1e; 1570 } else if (compare < 0) { 1571 #if BOOST_UBLAS_TYPE_CHECK 1572 while (it1.index1 () < it1e.index1 ()) { 1573 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1574 typename M::iterator2 it2 (it1.begin ()); 1575 typename M::iterator2 it2_end (it1.end ()); 1576 #else 1577 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 1578 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 1579 #endif 1580 increment (it2, it2_end); 1581 ++ it1; 1582 } 1583 #else 1584 increment (it1, it1_end, - compare); 1585 #endif 1586 } else if (compare > 0) { 1587 #if BOOST_UBLAS_TYPE_CHECK 1588 while (it1e.index1 () < it1.index1 ()) { 1589 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1590 typename E::iterator2 it2e (it1e.begin ()); 1591 typename E::iterator2 it2e_end (it1e.end ()); 1592 #else 1593 typename E::iterator2 it2e (begin (it1e, iterator1_tag ())); 1594 typename E::iterator2 it2e_end (end (it1e, iterator1_tag ())); 1595 #endif 1596 increment (it2e, it2e_end); 1597 ++ it1e; 1598 } 1599 #else 1600 increment (it1e, it1e_end, compare); 1601 #endif 1602 } 1603 } 1604 #if BOOST_UBLAS_TYPE_CHECK 1605 while (it1e != it1e_end) { 1606 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1607 typename E::iterator2 it2e (it1e.begin ()); 1608 typename E::iterator2 it2e_end (it1e.end ()); 1609 #else 1610 typename E::iterator2 it2e (begin (it1e, iterator1_tag ())); 1611 typename E::iterator2 it2e_end (end (it1e, iterator1_tag ())); 1612 #endif 1613 increment (it2e, it2e_end); 1614 ++ it1e; 1615 } 1616 while (it1 != it1_end) { 1617 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1618 typename M::iterator2 it2 (it1.begin ()); 1619 typename M::iterator2 it2_end (it1.end ()); 1620 #else 1621 typename M::iterator2 it2 (begin (it1, iterator1_tag ())); 1622 typename M::iterator2 it2_end (end (it1, iterator1_tag ())); 1623 #endif 1624 increment (it2, it2_end); 1625 ++ it1; 1626 } 1627 #endif 1628 } 1629 // Sparse (proxy) column major case 1630 template<template <class T1, class T2> class F, class R, class M, class E> 1631 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. matrix_swap(M & m,matrix_expression<E> & e,sparse_proxy_tag,column_major_tag)1632 void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) { 1633 typedef F<typename M::iterator1::reference, typename E::reference> functor_type; 1634 typedef R conformant_restrict_type; 1635 typedef typename M::size_type size_type; 1636 typedef typename M::difference_type difference_type; 1637 1638 BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ()); 1639 BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ()); 1640 1641 detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ()); 1642 // FIXME should be a seperate restriction for E 1643 detail::make_conformant (e (), m, column_major_tag (), conformant_restrict_type ()); 1644 1645 typename M::iterator2 it2 (m.begin2 ()); 1646 typename M::iterator2 it2_end (m.end2 ()); 1647 typename E::iterator2 it2e (e ().begin2 ()); 1648 typename E::iterator2 it2e_end (e ().end2 ()); 1649 while (it2 != it2_end && it2e != it2e_end) { 1650 difference_type compare = it2.index2 () - it2e.index2 (); 1651 if (compare == 0) { 1652 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1653 typename M::iterator1 it1 (it2.begin ()); 1654 typename M::iterator1 it1_end (it2.end ()); 1655 typename E::iterator1 it1e (it2e.begin ()); 1656 typename E::iterator1 it1e_end (it2e.end ()); 1657 #else 1658 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 1659 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 1660 typename E::iterator1 it1e (begin (it2e, iterator2_tag ())); 1661 typename E::iterator1 it1e_end (end (it2e, iterator2_tag ())); 1662 #endif 1663 if (it1 != it1_end && it1e != it1e_end) { 1664 size_type it1_index = it1.index1 (), it1e_index = it1e.index1 (); 1665 for (;;) { 1666 difference_type compare2 = it1_index - it1e_index; 1667 if (compare2 == 0) { 1668 functor_type::apply (*it1, *it1e); 1669 ++ it1, ++ it1e; 1670 if (it1 != it1_end && it1e != it1e_end) { 1671 it1_index = it1.index1 (); 1672 it1e_index = it1e.index1 (); 1673 } else 1674 break; 1675 } else if (compare2 < 0) { 1676 increment (it1, it1_end, - compare2); 1677 if (it1 != it1_end) 1678 it1_index = it1.index1 (); 1679 else 1680 break; 1681 } else if (compare2 > 0) { 1682 increment (it1e, it1e_end, compare2); 1683 if (it1e != it1e_end) 1684 it1e_index = it1e.index1 (); 1685 else 1686 break; 1687 } 1688 } 1689 } 1690 #if BOOST_UBLAS_TYPE_CHECK 1691 increment (it1e, it1e_end); 1692 increment (it1, it1_end); 1693 #endif 1694 ++ it2, ++ it2e; 1695 } else if (compare < 0) { 1696 #if BOOST_UBLAS_TYPE_CHECK 1697 while (it2.index2 () < it2e.index2 ()) { 1698 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1699 typename M::iterator1 it1 (it2.begin ()); 1700 typename M::iterator1 it1_end (it2.end ()); 1701 #else 1702 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 1703 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 1704 #endif 1705 increment (it1, it1_end); 1706 ++ it2; 1707 } 1708 #else 1709 increment (it2, it2_end, - compare); 1710 #endif 1711 } else if (compare > 0) { 1712 #if BOOST_UBLAS_TYPE_CHECK 1713 while (it2e.index2 () < it2.index2 ()) { 1714 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1715 typename E::iterator1 it1e (it2e.begin ()); 1716 typename E::iterator1 it1e_end (it2e.end ()); 1717 #else 1718 typename E::iterator1 it1e (begin (it2e, iterator2_tag ())); 1719 typename E::iterator1 it1e_end (end (it2e, iterator2_tag ())); 1720 #endif 1721 increment (it1e, it1e_end); 1722 ++ it2e; 1723 } 1724 #else 1725 increment (it2e, it2e_end, compare); 1726 #endif 1727 } 1728 } 1729 #if BOOST_UBLAS_TYPE_CHECK 1730 while (it2e != it2e_end) { 1731 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1732 typename E::iterator1 it1e (it2e.begin ()); 1733 typename E::iterator1 it1e_end (it2e.end ()); 1734 #else 1735 typename E::iterator1 it1e (begin (it2e, iterator2_tag ())); 1736 typename E::iterator1 it1e_end (end (it2e, iterator2_tag ())); 1737 #endif 1738 increment (it1e, it1e_end); 1739 ++ it2e; 1740 } 1741 while (it2 != it2_end) { 1742 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1743 typename M::iterator1 it1 (it2.begin ()); 1744 typename M::iterator1 it1_end (it2.end ()); 1745 #else 1746 typename M::iterator1 it1 (begin (it2, iterator2_tag ())); 1747 typename M::iterator1 it1_end (end (it2, iterator2_tag ())); 1748 #endif 1749 increment (it1, it1_end); 1750 ++ it2; 1751 } 1752 #endif 1753 } 1754 1755 // Dispatcher 1756 template<template <class T1, class T2> class F, class M, class E> 1757 BOOST_UBLAS_INLINE matrix_swap(M & m,matrix_expression<E> & e)1758 void matrix_swap (M &m, matrix_expression<E> &e) { 1759 typedef typename matrix_swap_traits<typename M::storage_category, 1760 typename E::const_iterator1::iterator_category, 1761 typename E::const_iterator2::iterator_category>::storage_category storage_category; 1762 // give preference to matrix M's orientation if known 1763 typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>, 1764 typename E::orientation_category , 1765 typename M::orientation_category >::type orientation_category; 1766 typedef basic_full<typename M::size_type> unrestricted; 1767 matrix_swap<F, unrestricted> (m, e, storage_category (), orientation_category ()); 1768 } 1769 template<template <class T1, class T2> class F, class R, class M, class E> 1770 BOOST_UBLAS_INLINE matrix_swap(M & m,matrix_expression<E> & e)1771 void matrix_swap (M &m, matrix_expression<E> &e) { 1772 typedef R conformant_restrict_type; 1773 typedef typename matrix_swap_traits<typename M::storage_category, 1774 typename E::const_iterator1::iterator_category, 1775 typename E::const_iterator2::iterator_category>::storage_category storage_category; 1776 // give preference to matrix M's orientation if known 1777 typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>, 1778 typename E::orientation_category , 1779 typename M::orientation_category >::type orientation_category; 1780 matrix_swap<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ()); 1781 } 1782 1783 }}} 1784 1785 #endif 1786