boost/numeric/ublas/detail/matrix_assign.hpp

   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
  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
  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.
  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                     while (true) {
  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.
 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                     while (true) {
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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
 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.
 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.
 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.
 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.
 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.
 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.
 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                 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 ());
 678 #else
 679                         typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
 680                         typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
 681 #endif
 682                         difference_type size2 (it2_end - it2);
 683                         while (-- size2 >= 0)
 684                             functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
 685                         ++ it1;
 686                     }
 687                 } else {
 688                     it1 += size1;
 689                 }
 690                 diff1 += size1;
 691             }
 692         }
 693         difference_type size1 ((std::min) (it1_size, it1e_size));
 694         it1_size -= size1;
 695         it1e_size -= size1;
 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 ());
 702 #else
 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 ()));
 707 #endif
 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);
 714                 if (size2 > 0) {
 715                     it2e += size2;
 716                     it2e_size -= size2;
 717                     diff2 -= size2;
 718                 }
 719                 size2 = (std::min) (- diff2, it2_size);
 720                 if (size2 > 0) {
 721                     it2_size -= size2;
 722                     if (!functor_type::computed) {
 723                         while (-- size2 >= 0)   // zeroing
 724                             functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
 725                     } else {
 726                         it2 += size2;
 727                     }
 728                     diff2 += size2;
 729                 }
 730             }
 731             difference_type size2 ((std::min) (it2_size, it2e_size));
 732             it2_size -= size2;
 733             it2e_size -= size2;
 734             while (-- size2 >= 0)
 735                 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
 736             size2 = it2_size;
 737             if (!functor_type::computed) {
 738                 while (-- size2 >= 0)   // zeroing
 739                     functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
 740             } else {
 741                 it2 += size2;
 742             }
 743             ++ it1, ++ it1e;
 744         }
 745         size1 = it1_size;
 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 ());
 751 #else
 752                 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
 753                 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
 754 #endif
 755                 difference_type size2 (it2_end - it2);
 756                 while (-- size2 >= 0)
 757                     functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
 758                 ++ it1;
 759             }
 760         } else {
 761             it1 += size1;
 762         }
 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 ());
 766 #endif
 767     }
 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;
 776
 777         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
 778         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
 779
 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 ());
 785 #endif
 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 ();
 795         if (diff2 != 0) {
 796             difference_type size2 = (std::min) (diff2, it2e_size);
 797             if (size2 > 0) {
 798                 it2e += size2;
 799                 it2e_size -= size2;
 800                 diff2 -= size2;
 801             }
 802             size2 = (std::min) (- diff2, it2_size);
 803             if (size2 > 0) {
 804                 it2_size -= size2;
 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 ());
 810 #else
 811                         typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
 812                         typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
 813 #endif
 814                         difference_type size1 (it1_end - it1);
 815                         while (-- size1 >= 0)
 816                             functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
 817                         ++ it2;
 818                     }
 819                 } else {
 820                     it2 += size2;
 821                 }
 822                 diff2 += size2;
 823             }
 824         }
 825         difference_type size2 ((std::min) (it2_size, it2e_size));
 826         it2_size -= size2;
 827         it2e_size -= size2;
 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 ());
 834 #else
 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 ()));
 839 #endif
 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);
 846                 if (size1 > 0) {
 847                     it1e += size1;
 848                     it1e_size -= size1;
 849                     diff1 -= size1;
 850                 }
 851                 size1 = (std::min) (- diff1, it1_size);
 852                 if (size1 > 0) {
 853                     it1_size -= size1;
 854                     if (!functor_type::computed) {
 855                         while (-- size1 >= 0)   // zeroing
 856                             functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
 857                     } else {
 858                         it1 += size1;
 859                     }
 860                     diff1 += size1;
 861                 }
 862             }
 863             difference_type size1 ((std::min) (it1_size, it1e_size));
 864             it1_size -= size1;
 865             it1e_size -= size1;
 866             while (-- size1 >= 0)
 867                 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
 868             size1 = it1_size;
 869             if (!functor_type::computed) {
 870                 while (-- size1 >= 0)   // zeroing
 871                     functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
 872             } else {
 873                 it1 += size1;
 874             }
 875             ++ it2, ++ it2e;
 876         }
 877         size2 = it2_size;
 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 ());
 883 #else
 884                 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
 885                 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
 886 #endif
 887                 difference_type size1 (it1_end - it1);
 888                 while (-- size1 >= 0)
 889                     functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
 890                 ++ it2;
 891             }
 892         } else {
 893             it2 += size2;
 894         }
 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 ());
 898 #endif
 899     }
 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
 911
 912         m.clear ();
 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 ());
 919 #else
 920             typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
 921             typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
 922 #endif
 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);
 927                 ++ it2e;
 928             }
 929             ++ it1e;
 930         }
 931     }
 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
 943
 944         m.clear ();
 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 ());
 951 #else
 952             typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
 953             typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
 954 #endif
 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);
 959                 ++ it1e;
 960             }
 961             ++ it2e;
 962         }
 963     }
 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;
 973
 974         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
 975         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
 976
 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 ());
 982 #endif
 983         detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
 984
 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 ();
 991             if (compare == 0) {
 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 ());
 997 #else
 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 ()));
1002 #endif
1003                 if (it2 != it2_end && it2e != it2e_end) {
1004                     size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
1005                     while (true) {
1006                         difference_type compare2 = it2_index - it2e_index;
1007                         if (compare2 == 0) {
1008                             functor_type::apply (*it2, *it2e);
1009                             ++ it2, ++ it2e;
1010                             if (it2 != it2_end && it2e != it2e_end) {
1011                                 it2_index = it2.index2 ();
1012                                 it2e_index = it2e.index2 ();
1013                             } else
1014                                 break;
1015                         } else if (compare2 < 0) {
1016                             if (!functor_type::computed) {
1017                                 functor_type::apply (*it2, expr_value_type/*zero*/());
1018                                 ++ it2;
1019                             } else
1020                                 increment (it2, it2_end, - compare2);
1021                             if (it2 != it2_end)
1022                                 it2_index = it2.index2 ();
1023                             else
1024                                 break;
1025                         } else if (compare2 > 0) {
1026                             increment (it2e, it2e_end, compare2);
1027                             if (it2e != it2e_end)
1028                                 it2e_index = it2e.index2 ();
1029                             else
1030                                 break;
1031                         }
1032                     }
1033                 }
1034                 if (!functor_type::computed) {
1035                     while (it2 != it2_end) {    // zeroing
1036                         functor_type::apply (*it2, expr_value_type/*zero*/());
1037                         ++ it2;
1038                     }
1039                 } else {
1040                     it2 = it2_end;
1041                 }
1042                 ++ it1, ++ it1e;
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 ());
1048 #else
1049                     typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1050                     typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1051 #endif
1052                     while (it2 != it2_end) {    // zeroing
1053                         functor_type::apply (*it2, expr_value_type/*zero*/());
1054                         ++ it2;
1055                     }
1056                     ++ it1;
1057                 } else {
1058                     increment (it1, it1_end, - compare);
1059                 }
1060             } else if (compare > 0) {
1061                 increment (it1e, it1e_end, compare);
1062             }
1063         }
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 ());
1069 #else
1070                 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1071                 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1072 #endif
1073                 while (it2 != it2_end) {    // zeroing
1074                     functor_type::apply (*it2, expr_value_type/*zero*/());
1075                     ++ it2;
1076                 }
1077                 ++ it1;
1078             }
1079         } else {
1080             it1 = it1_end;
1081         }
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 ());
1085 #endif
1086     }
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;
1096
1097         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
1098         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
1099
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 ());
1105 #endif
1106         detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
1107
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 ();
1114             if (compare == 0) {
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 ());
1120 #else
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 ()));
1125 #endif
1126                 if (it1 != it1_end && it1e != it1e_end) {
1127                     size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
1128                     while (true) {
1129                         difference_type compare2 = it1_index - it1e_index;
1130                         if (compare2 == 0) {
1131                             functor_type::apply (*it1, *it1e);
1132                             ++ it1, ++ it1e;
1133                             if (it1 != it1_end && it1e != it1e_end) {
1134                                 it1_index = it1.index1 ();
1135                                 it1e_index = it1e.index1 ();
1136                             } else
1137                                 break;
1138                         } else if (compare2 < 0) {
1139                             if (!functor_type::computed) {
1140                                 functor_type::apply (*it1, expr_value_type/*zero*/()); // zeroing
1141                                 ++ it1;
1142                             } else
1143                                 increment (it1, it1_end, - compare2);
1144                             if (it1 != it1_end)
1145                                 it1_index = it1.index1 ();
1146                             else
1147                                 break;
1148                         } else if (compare2 > 0) {
1149                             increment (it1e, it1e_end, compare2);
1150                             if (it1e != it1e_end)
1151                                 it1e_index = it1e.index1 ();
1152                             else
1153                                 break;
1154                         }
1155                     }
1156                 }
1157                 if (!functor_type::computed) {
1158                     while (it1 != it1_end) {    // zeroing
1159                         functor_type::apply (*it1, expr_value_type/*zero*/());
1160                         ++ it1;
1161                     }
1162                 } else {
1163                     it1 = it1_end;
1164                 }
1165                 ++ it2, ++ it2e;
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 ());
1171 #else
1172                     typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1173                     typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1174 #endif
1175                     while (it1 != it1_end) {    // zeroing
1176                         functor_type::apply (*it1, expr_value_type/*zero*/());
1177                         ++ it1;
1178                     }
1179                     ++ it2;
1180                 } else {
1181                     increment (it2, it2_end, - compare);
1182                 }
1183             } else if (compare > 0) {
1184                 increment (it2e, it2e_end, compare);
1185             }
1186         }
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 ());
1192 #else
1193                 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1194                 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1195 #endif
1196                 while (it1 != it1_end) {    // zeroing
1197                     functor_type::apply (*it1, expr_value_type/*zero*/());
1198                     ++ it1;
1199                 }
1200                 ++ it2;
1201             }
1202         } else {
1203             it2 = it2_end;
1204         }
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 ());
1208 #endif
1209     }
1210
1211     // Dispatcher
1212     template<template <class T1, class T2> class F, class M, class E>
1213     BOOST_UBLAS_INLINE
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 ());
1225     }
1226     template<template <class T1, class T2> class F, class R, class M, class E>
1227     BOOST_UBLAS_INLINE
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 ());
1239     }
1240
1241     template<class SC, class RI1, class RI2>
1242     struct matrix_swap_traits {
1243         typedef SC storage_category;
1244     };
1245
1246     template<>
1247     struct matrix_swap_traits<dense_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
1248         typedef sparse_proxy_tag storage_category;
1249     };
1250
1251     template<>
1252     struct matrix_swap_traits<packed_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
1253         typedef sparse_proxy_tag storage_category;
1254     };
1255
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;
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 (), 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 (), size_type (it1e.end () - it2e)));
1272 #else
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 (), size_type (end (it1e, iterator1_tag ()) - it2e)));
1276 #endif
1277             while (-- size2 >= 0)
1278                 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
1279             ++ it1, ++ it1e;
1280         }
1281     }
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;
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 (), 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 (), size_type (it2e.end () - it1e)));
1298 #else
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 (), size_type (end (it2e, iterator2_tag ()) - it1e)));
1302 #endif
1303             while (-- size1 >= 0)
1304                 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
1305             ++ it2, ++ it2e;
1306         }
1307     }
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));
1323 #else
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));
1327 #endif
1328             while (-- size2 >= 0)
1329                 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
1330             ++ it1, ++ it1e;
1331         }
1332     }
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));
1348 #else
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));
1352 #endif
1353             while (-- size1 >= 0)
1354                 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
1355             ++ it2, ++ it2e;
1356         }
1357     }
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 ());
1368
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 ());
1372
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 ();
1379             if (compare == 0) {
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 ());
1385 #else
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 ()));
1390 #endif
1391                 if (it2 != it2_end && it2e != it2e_end) {
1392                     size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
1393                     while (true) {
1394                         difference_type compare2 = it2_index - it2e_index;
1395                         if (compare2 == 0) {
1396                             functor_type::apply (*it2, *it2e);
1397                             ++ it2, ++ it2e;
1398                             if (it2 != it2_end && it2e != it2e_end) {
1399                                 it2_index = it2.index2 ();
1400                                 it2e_index = it2e.index2 ();
1401                             } else
1402                                 break;
1403                         } else if (compare2 < 0) {
1404                             increment (it2, it2_end, - compare2);
1405                             if (it2 != it2_end)
1406                                 it2_index = it2.index2 ();
1407                             else
1408                                 break;
1409                         } else if (compare2 > 0) {
1410                             increment (it2e, it2e_end, compare2);
1411                             if (it2e != it2e_end)
1412                                 it2e_index = it2e.index2 ();
1413                             else
1414                                 break;
1415                         }
1416                     }
1417                 }
1418 #if BOOST_UBLAS_TYPE_CHECK
1419                 increment (it2e, it2e_end);
1420                 increment (it2, it2_end);
1421 #endif
1422                 ++ it1, ++ it1e;
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 ());
1429 #else
1430                     typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1431                     typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1432 #endif
1433                     increment (it2, it2_end);
1434                     ++ it1;
1435                 }
1436 #else
1437                 increment (it1, it1_end, - compare);
1438 #endif
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 ());
1445 #else
1446                     typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1447                     typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
1448 #endif
1449                     increment (it2e, it2e_end);
1450                     ++ it1e;
1451                 }
1452 #else
1453                 increment (it1e, it1e_end, compare);
1454 #endif
1455             }
1456         }
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 ());
1462 #else
1463             typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1464             typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
1465 #endif
1466             increment (it2e, it2e_end);
1467             ++ it1e;
1468         }
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 ());
1473 #else
1474             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1475             typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1476 #endif
1477             increment (it2, it2_end);
1478             ++ it1;
1479         }
1480 #endif
1481     }
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;
1490
1491         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
1492         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
1493
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 ());
1497
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 ();
1504             if (compare == 0) {
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 ());
1510 #else
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 ()));
1515 #endif
1516                 if (it1 != it1_end && it1e != it1e_end) {
1517                     size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
1518                     while (true) {
1519                         difference_type compare2 = it1_index - it1e_index;
1520                         if (compare2 == 0) {
1521                             functor_type::apply (*it1, *it1e);
1522                             ++ it1, ++ it1e;
1523                             if (it1 != it1_end && it1e != it1e_end) {
1524                                 it1_index = it1.index1 ();
1525                                 it1e_index = it1e.index1 ();
1526                             } else
1527                                 break;
1528                         }  else if (compare2 < 0) {
1529                             increment (it1, it1_end, - compare2);
1530                             if (it1 != it1_end)
1531                                 it1_index = it1.index1 ();
1532                             else
1533                                 break;
1534                         } else if (compare2 > 0) {
1535                             increment (it1e, it1e_end, compare2);
1536                             if (it1e != it1e_end)
1537                                 it1e_index = it1e.index1 ();
1538                             else
1539                                 break;
1540                         }
1541                     }
1542                 }
1543 #if BOOST_UBLAS_TYPE_CHECK
1544                 increment (it1e, it1e_end);
1545                 increment (it1, it1_end);
1546 #endif
1547                 ++ it2, ++ it2e;
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 ());
1554 #else
1555                     typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1556                     typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1557 #endif
1558                     increment (it1, it1_end);
1559                     ++ it2;
1560                 }
1561 #else
1562                 increment (it2, it2_end, - compare);
1563 #endif
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 ());
1570 #else
1571                     typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1572                     typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
1573 #endif
1574                     increment (it1e, it1e_end);
1575                     ++ it2e;
1576                 }
1577 #else
1578                 increment (it2e, it2e_end, compare);
1579 #endif
1580             }
1581         }
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 ());
1587 #else
1588             typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1589             typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
1590 #endif
1591             increment (it1e, it1e_end);
1592             ++ it2e;
1593         }
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 ());
1598 #else
1599             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1600             typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1601 #endif
1602             increment (it1, it1_end);
1603             ++ it2;
1604         }
1605 #endif
1606     }
1607
1608     // Dispatcher
1609     template<template <class T1, class T2> class F, class M, class E>
1610     BOOST_UBLAS_INLINE
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 ());
1621     }
1622     template<template <class T1, class T2> class F, class R, class M, class E>
1623     BOOST_UBLAS_INLINE
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 ());
1634     }
1635
1636 }}}
1637
1638 #endif