ceph/src/boost/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                     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.
 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.
 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 //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.
 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.
 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.
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.
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.
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
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
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.
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.
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.
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.
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.
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.
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
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
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