Boost.uBlas 1.49
Linear Algebra in C++: matrices, vectors and numeric algorithms

matrix_assign.hpp

Go to the documentation of this file.
00001 //
00002 //  Copyright (c) 2000-2002
00003 //  Joerg Walter, Mathias Koch
00004 //
00005 //  Distributed under the Boost Software License, Version 1.0. (See
00006 //  accompanying file LICENSE_1_0.txt or copy at
00007 //  http://www.boost.org/LICENSE_1_0.txt)
00008 //
00009 //  The authors gratefully acknowledge the support of
00010 //  GeNeSys mbH & Co. KG in producing this work.
00011 //
00012 
00013 #ifndef _BOOST_UBLAS_MATRIX_ASSIGN_
00014 #define _BOOST_UBLAS_MATRIX_ASSIGN_
00015 
00016 #include <boost/numeric/ublas/traits.hpp>
00017 // Required for make_conformant storage
00018 #include <vector>
00019 
00020 // Iterators based on ideas of Jeremy Siek
00021 
00022 namespace boost { namespace numeric { namespace ublas {
00023 namespace detail {
00024     
00025     // Weak equality check - useful to compare equality two arbitary matrix expression results.
00026     // Since the actual expressions are unknown, we check for and arbitary error bound
00027     // on the relative error.
00028     // For a linear expression the infinity norm makes sense as we do not know how the elements will be
00029     // combined in the expression. False positive results are inevitable for arbirary expressions!
00030     template<class E1, class E2, class S>
00031     BOOST_UBLAS_INLINE
00032     bool equals (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2, S epsilon, S min_norm) {
00033         return norm_inf (e1 - e2) < epsilon *
00034                std::max<S> (std::max<S> (norm_inf (e1), norm_inf (e2)), min_norm);
00035     }
00036 
00037     template<class E1, class E2>
00038     BOOST_UBLAS_INLINE
00039     bool expression_type_check (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2) {
00040         typedef typename type_traits<typename promote_traits<typename E1::value_type,
00041                                      typename E2::value_type>::promote_type>::real_type real_type;
00042         return equals (e1, e2, BOOST_UBLAS_TYPE_CHECK_EPSILON, BOOST_UBLAS_TYPE_CHECK_MIN);
00043     }
00044 
00045 
00046     template<class M, class E, class R>
00047     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00048     void make_conformant (M &m, const matrix_expression<E> &e, row_major_tag, R) {
00049         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
00050         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
00051         typedef R conformant_restrict_type;
00052         typedef typename M::size_type size_type;
00053         typedef typename M::difference_type difference_type;
00054         typedef typename M::value_type value_type;
00055         // FIXME unbounded_array with push_back maybe better
00056         std::vector<std::pair<size_type, size_type> > index;
00057         typename M::iterator1 it1 (m.begin1 ());
00058         typename M::iterator1 it1_end (m.end1 ());
00059         typename E::const_iterator1 it1e (e ().begin1 ());
00060         typename E::const_iterator1 it1e_end (e ().end1 ());
00061         while (it1 != it1_end && it1e != it1e_end) {
00062             difference_type compare = it1.index1 () - it1e.index1 ();
00063             if (compare == 0) {
00064 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00065                 typename M::iterator2 it2 (it1.begin ());
00066                 typename M::iterator2 it2_end (it1.end ());
00067                 typename E::const_iterator2 it2e (it1e.begin ());
00068                 typename E::const_iterator2 it2e_end (it1e.end ());
00069 #else
00070                 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00071                 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
00072                 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
00073                 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
00074 #endif
00075                 if (it2 != it2_end && it2e != it2e_end) {
00076                     size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
00077                     while (true) {
00078                         difference_type compare = it2_index - it2e_index;
00079                         if (compare == 0) {
00080                             ++ it2, ++ it2e;
00081                             if (it2 != it2_end && it2e != it2e_end) {
00082                                 it2_index = it2.index2 ();
00083                                 it2e_index = it2e.index2 ();
00084                             } else
00085                                 break;
00086                         } else if (compare < 0) {
00087                             increment (it2, it2_end, - compare);
00088                             if (it2 != it2_end)
00089                                 it2_index = it2.index2 ();
00090                             else
00091                                 break;
00092                         } else if (compare > 0) {
00093                             if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
00094                                 if (static_cast<value_type>(*it2e) != value_type/*zero*/())
00095                                     index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
00096                             ++ it2e;
00097                             if (it2e != it2e_end)
00098                                 it2e_index = it2e.index2 ();
00099                             else
00100                                 break;
00101                         }
00102                     }
00103                 }
00104                 while (it2e != it2e_end) {
00105                     if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
00106                         if (static_cast<value_type>(*it2e) != value_type/*zero*/())
00107                             index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
00108                     ++ it2e;
00109                 }
00110                 ++ it1, ++ it1e;
00111             } else if (compare < 0) {
00112                 increment (it1, it1_end, - compare);
00113             } else if (compare > 0) {
00114 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00115                 typename E::const_iterator2 it2e (it1e.begin ());
00116                 typename E::const_iterator2 it2e_end (it1e.end ());
00117 #else
00118                 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
00119                 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
00120 #endif
00121                 while (it2e != it2e_end) {
00122                     if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
00123                         if (static_cast<value_type>(*it2e) != value_type/*zero*/())
00124                             index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
00125                     ++ it2e;
00126                 }
00127                 ++ it1e;
00128             }
00129         }
00130         while (it1e != it1e_end) {
00131 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00132             typename E::const_iterator2 it2e (it1e.begin ());
00133             typename E::const_iterator2 it2e_end (it1e.end ());
00134 #else
00135             typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
00136             typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
00137 #endif
00138             while (it2e != it2e_end) {
00139                 if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
00140                     if (static_cast<value_type>(*it2e) != value_type/*zero*/())
00141                         index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
00142                 ++ it2e;
00143             }
00144             ++ it1e;
00145         }
00146         // ISSUE proxies require insert_element
00147         for (size_type k = 0; k < index.size (); ++ k)
00148             m (index [k].first, index [k].second) = value_type/*zero*/();
00149     }
00150     template<class M, class E, class R>
00151     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00152     void make_conformant (M &m, const matrix_expression<E> &e, column_major_tag, R) {
00153         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
00154         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
00155         typedef R conformant_restrict_type;
00156         typedef typename M::size_type size_type;
00157         typedef typename M::difference_type difference_type;
00158         typedef typename M::value_type value_type;
00159         std::vector<std::pair<size_type, size_type> > index;
00160         typename M::iterator2 it2 (m.begin2 ());
00161         typename M::iterator2 it2_end (m.end2 ());
00162         typename E::const_iterator2 it2e (e ().begin2 ());
00163         typename E::const_iterator2 it2e_end (e ().end2 ());
00164         while (it2 != it2_end && it2e != it2e_end) {
00165             difference_type compare = it2.index2 () - it2e.index2 ();
00166             if (compare == 0) {
00167 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00168                 typename M::iterator1 it1 (it2.begin ());
00169                 typename M::iterator1 it1_end (it2.end ());
00170                 typename E::const_iterator1 it1e (it2e.begin ());
00171                 typename E::const_iterator1 it1e_end (it2e.end ());
00172 #else
00173                 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
00174                 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
00175                 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
00176                 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
00177 #endif
00178                 if (it1 != it1_end && it1e != it1e_end) {
00179                     size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
00180                     while (true) {
00181                         difference_type compare = it1_index - it1e_index;
00182                         if (compare == 0) {
00183                             ++ it1, ++ it1e;
00184                             if (it1 != it1_end && it1e != it1e_end) {
00185                                 it1_index = it1.index1 ();
00186                                 it1e_index = it1e.index1 ();
00187                             } else
00188                                 break;
00189                         } else if (compare < 0) {
00190                             increment (it1, it1_end, - compare);
00191                             if (it1 != it1_end)
00192                                 it1_index = it1.index1 ();
00193                             else
00194                                 break;
00195                         } else if (compare > 0) {
00196                             if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
00197                                 if (static_cast<value_type>(*it1e) != value_type/*zero*/())
00198                                     index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
00199                             ++ it1e;
00200                             if (it1e != it1e_end)
00201                                 it1e_index = it1e.index1 ();
00202                             else
00203                                 break;
00204                         }
00205                     }
00206                 }
00207                 while (it1e != it1e_end) {
00208                     if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
00209                         if (static_cast<value_type>(*it1e) != value_type/*zero*/())
00210                             index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
00211                     ++ it1e;
00212                 }
00213                 ++ it2, ++ it2e;
00214             } else if (compare < 0) {
00215                 increment (it2, it2_end, - compare);
00216             } else if (compare > 0) {
00217 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00218                 typename E::const_iterator1 it1e (it2e.begin ());
00219                 typename E::const_iterator1 it1e_end (it2e.end ());
00220 #else
00221                 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
00222                 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
00223 #endif
00224                 while (it1e != it1e_end) {
00225                     if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
00226                         if (static_cast<value_type>(*it1e) != value_type/*zero*/())
00227                             index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
00228                     ++ it1e;
00229                 }
00230                 ++ it2e;
00231             }
00232         }
00233         while (it2e != it2e_end) {
00234 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00235             typename E::const_iterator1 it1e (it2e.begin ());
00236             typename E::const_iterator1 it1e_end (it2e.end ());
00237 #else
00238             typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
00239             typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
00240 #endif
00241             while (it1e != it1e_end) {
00242                 if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
00243                     if (static_cast<value_type>(*it1e) != value_type/*zero*/())
00244                         index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
00245                 ++ it1e;
00246             }
00247             ++ it2e;
00248         }
00249         // ISSUE proxies require insert_element
00250         for (size_type k = 0; k < index.size (); ++ k)
00251             m (index [k].first, index [k].second) = value_type/*zero*/();
00252     }
00253 
00254 }//namespace detail
00255 
00256 
00257     // Explicitly iterating row major
00258     template<template <class T1, class T2> class F, class M, class T>
00259     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00260     void iterating_matrix_assign_scalar (M &m, const T &t, row_major_tag) {
00261         typedef F<typename M::iterator2::reference, T> functor_type;
00262         typedef typename M::difference_type difference_type;
00263         difference_type size1 (m.size1 ());
00264         difference_type size2 (m.size2 ());
00265         typename M::iterator1 it1 (m.begin1 ());
00266         BOOST_UBLAS_CHECK (size2 == 0 || m.end1 () - it1 == size1, bad_size ());
00267         while (-- size1 >= 0) {
00268 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00269             typename M::iterator2 it2 (it1.begin ());
00270 #else
00271             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00272 #endif
00273             BOOST_UBLAS_CHECK (it1.end () - it2 == size2, bad_size ());
00274             difference_type temp_size2 (size2);
00275 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
00276             while (-- temp_size2 >= 0)
00277                 functor_type::apply (*it2, t), ++ it2;
00278 #else
00279             DD (temp_size2, 4, r, (functor_type::apply (*it2, t), ++ it2));
00280 #endif
00281             ++ it1;
00282         }
00283     }
00284     // Explicitly iterating column major
00285     template<template <class T1, class T2> class F, class M, class T>
00286     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00287     void iterating_matrix_assign_scalar (M &m, const T &t, column_major_tag) {
00288         typedef F<typename M::iterator1::reference, T> functor_type;
00289         typedef typename M::difference_type difference_type;
00290         difference_type size2 (m.size2 ());
00291         difference_type size1 (m.size1 ());
00292         typename M::iterator2 it2 (m.begin2 ());
00293         BOOST_UBLAS_CHECK (size1 == 0 || m.end2 () - it2 == size2, bad_size ());
00294         while (-- size2 >= 0) {
00295 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00296             typename M::iterator1 it1 (it2.begin ());
00297 #else
00298             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
00299 #endif
00300             BOOST_UBLAS_CHECK (it2.end () - it1 == size1, bad_size ());
00301             difference_type temp_size1 (size1);
00302 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
00303             while (-- temp_size1 >= 0)
00304                 functor_type::apply (*it1, t), ++ it1;
00305 #else
00306             DD (temp_size1, 4, r, (functor_type::apply (*it1, t), ++ it1));
00307 #endif
00308             ++ it2;
00309         }
00310     }
00311     // Explicitly indexing row major
00312     template<template <class T1, class T2> class F, class M, class T>
00313     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00314     void indexing_matrix_assign_scalar (M &m, const T &t, row_major_tag) {
00315         typedef F<typename M::reference, T> functor_type;
00316         typedef typename M::size_type size_type;
00317         size_type size1 (m.size1 ());
00318         size_type size2 (m.size2 ());
00319         for (size_type i = 0; i < size1; ++ i) {
00320 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
00321             for (size_type j = 0; j < size2; ++ j)
00322                 functor_type::apply (m (i, j), t);
00323 #else
00324             size_type j (0);
00325             DD (size2, 4, r, (functor_type::apply (m (i, j), t), ++ j));
00326 #endif
00327         }
00328     }
00329     // Explicitly indexing column major
00330     template<template <class T1, class T2> class F, class M, class T>
00331     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00332     void indexing_matrix_assign_scalar (M &m, const T &t, column_major_tag) {
00333         typedef F<typename M::reference, T> functor_type;
00334         typedef typename M::size_type size_type;
00335         size_type size2 (m.size2 ());
00336         size_type size1 (m.size1 ());
00337         for (size_type j = 0; j < size2; ++ j) {
00338 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
00339             for (size_type i = 0; i < size1; ++ i)
00340                 functor_type::apply (m (i, j), t);
00341 #else
00342             size_type i (0);
00343             DD (size1, 4, r, (functor_type::apply (m (i, j), t), ++ i));
00344 #endif
00345         }
00346     }
00347 
00348     // Dense (proxy) case
00349     template<template <class T1, class T2> class F, class M, class T, class C>
00350     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00351     void matrix_assign_scalar (M &m, const T &t, dense_proxy_tag, C) {
00352         typedef C orientation_category;
00353 #ifdef BOOST_UBLAS_USE_INDEXING
00354         indexing_matrix_assign_scalar<F> (m, t, orientation_category ());
00355 #elif BOOST_UBLAS_USE_ITERATING
00356         iterating_matrix_assign_scalar<F> (m, t, orientation_category ());
00357 #else
00358         typedef typename M::size_type size_type;
00359         size_type size1 (m.size1 ());
00360         size_type size2 (m.size2 ());
00361         if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
00362             size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
00363             iterating_matrix_assign_scalar<F> (m, t, orientation_category ());
00364         else
00365             indexing_matrix_assign_scalar<F> (m, t, orientation_category ());
00366 #endif
00367     }
00368     // Packed (proxy) row major case
00369     template<template <class T1, class T2> class F, class M, class T>
00370     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00371     void matrix_assign_scalar (M &m, const T &t, packed_proxy_tag, row_major_tag) {
00372         typedef F<typename M::iterator2::reference, T> functor_type;
00373         typedef typename M::difference_type difference_type;
00374         typename M::iterator1 it1 (m.begin1 ());
00375         difference_type size1 (m.end1 () - it1);
00376         while (-- size1 >= 0) {
00377 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00378             typename M::iterator2 it2 (it1.begin ());
00379             difference_type size2 (it1.end () - it2);
00380 #else
00381             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00382             difference_type size2 (end (it1, iterator1_tag ()) - it2);
00383 #endif
00384             while (-- size2 >= 0)
00385                 functor_type::apply (*it2, t), ++ it2;
00386             ++ it1;
00387         }
00388     }
00389     // Packed (proxy) column major case
00390     template<template <class T1, class T2> class F, class M, class T>
00391     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00392     void matrix_assign_scalar (M &m, const T &t, packed_proxy_tag, column_major_tag) {
00393         typedef F<typename M::iterator1::reference, T> functor_type;
00394         typedef typename M::difference_type difference_type;
00395         typename M::iterator2 it2 (m.begin2 ());
00396         difference_type size2 (m.end2 () - it2);
00397         while (-- size2 >= 0) {
00398 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00399             typename M::iterator1 it1 (it2.begin ());
00400             difference_type size1 (it2.end () - it1);
00401 #else
00402             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
00403             difference_type size1 (end (it2, iterator2_tag ()) - it1);
00404 #endif
00405             while (-- size1 >= 0)
00406                 functor_type::apply (*it1, t), ++ it1;
00407             ++ it2;
00408         }
00409     }
00410     // Sparse (proxy) row major case
00411     template<template <class T1, class T2> class F, class M, class T>
00412     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00413     void matrix_assign_scalar (M &m, const T &t, sparse_proxy_tag, row_major_tag) {
00414         typedef F<typename M::iterator2::reference, T> functor_type;
00415         typename M::iterator1 it1 (m.begin1 ());
00416         typename M::iterator1 it1_end (m.end1 ());
00417         while (it1 != it1_end) {
00418 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00419             typename M::iterator2 it2 (it1.begin ());
00420             typename M::iterator2 it2_end (it1.end ());
00421 #else
00422             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00423             typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
00424 #endif
00425             while (it2 != it2_end)
00426                 functor_type::apply (*it2, t), ++ it2;
00427             ++ it1;
00428         }
00429     }
00430     // Sparse (proxy) column major case
00431     template<template <class T1, class T2> class F, class M, class T>
00432     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00433     void matrix_assign_scalar (M &m, const T &t, sparse_proxy_tag, column_major_tag) {
00434         typedef F<typename M::iterator1::reference, T> functor_type;
00435         typename M::iterator2 it2 (m.begin2 ());
00436         typename M::iterator2 it2_end (m.end2 ());
00437         while (it2 != it2_end) {
00438 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00439             typename M::iterator1 it1 (it2.begin ());
00440             typename M::iterator1 it1_end (it2.end ());
00441 #else
00442             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
00443             typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
00444 #endif
00445             while (it1 != it1_end)
00446                 functor_type::apply (*it1, t), ++ it1;
00447             ++ it2;
00448         }
00449     }
00450 
00451     // Dispatcher
00452     template<template <class T1, class T2> class F, class M, class T>
00453     BOOST_UBLAS_INLINE
00454     void matrix_assign_scalar (M &m, const T &t) {
00455         typedef typename M::storage_category storage_category;
00456         typedef typename M::orientation_category orientation_category;
00457         matrix_assign_scalar<F> (m, t, storage_category (), orientation_category ());
00458     }
00459 
00460     template<class SC, bool COMPUTED, class RI1, class RI2>
00461     struct matrix_assign_traits {
00462         typedef SC storage_category;
00463     };
00464 
00465     template<bool COMPUTED>
00466     struct matrix_assign_traits<dense_tag, COMPUTED, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
00467         typedef packed_tag storage_category;
00468     };
00469     template<>
00470     struct matrix_assign_traits<dense_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
00471         typedef sparse_tag storage_category;
00472     };
00473     template<>
00474     struct matrix_assign_traits<dense_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
00475         typedef sparse_proxy_tag storage_category;
00476     };
00477 
00478     template<bool COMPUTED>
00479     struct matrix_assign_traits<dense_proxy_tag, COMPUTED, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
00480         typedef packed_proxy_tag storage_category;
00481     };
00482     template<bool COMPUTED>
00483     struct matrix_assign_traits<dense_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
00484         typedef sparse_proxy_tag storage_category;
00485     };
00486 
00487     template<>
00488     struct matrix_assign_traits<packed_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
00489         typedef sparse_tag storage_category;
00490     };
00491     template<>
00492     struct matrix_assign_traits<packed_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
00493         typedef sparse_proxy_tag storage_category;
00494     };
00495 
00496     template<bool COMPUTED>
00497     struct matrix_assign_traits<packed_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
00498         typedef sparse_proxy_tag storage_category;
00499     };
00500 
00501     template<>
00502     struct matrix_assign_traits<sparse_tag, true, dense_random_access_iterator_tag, dense_random_access_iterator_tag> {
00503         typedef sparse_proxy_tag storage_category;
00504     };
00505     template<>
00506     struct matrix_assign_traits<sparse_tag, true, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
00507         typedef sparse_proxy_tag storage_category;
00508     };
00509     template<>
00510     struct matrix_assign_traits<sparse_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
00511         typedef sparse_proxy_tag storage_category;
00512     };
00513 
00514     // Explicitly iterating row major
00515     template<template <class T1, class T2> class F, class M, class E>
00516     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00517     void iterating_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
00518         typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
00519         typedef typename M::difference_type difference_type;
00520         difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
00521         difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
00522         typename M::iterator1 it1 (m.begin1 ());
00523         BOOST_UBLAS_CHECK (size2 == 0 || m.end1 () - it1 == size1, bad_size ());
00524         typename E::const_iterator1 it1e (e ().begin1 ());
00525         BOOST_UBLAS_CHECK (size2 == 0 || e ().end1 () - it1e == size1, bad_size ());
00526         while (-- size1 >= 0) {
00527 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00528             typename M::iterator2 it2 (it1.begin ());
00529             typename E::const_iterator2 it2e (it1e.begin ());
00530 #else
00531             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00532             typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
00533 #endif
00534             BOOST_UBLAS_CHECK (it1.end () - it2 == size2, bad_size ());
00535             BOOST_UBLAS_CHECK (it1e.end () - it2e == size2, bad_size ());
00536             difference_type temp_size2 (size2);
00537 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
00538             while (-- temp_size2 >= 0)
00539                 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
00540 #else
00541             DD (temp_size2, 2, r, (functor_type::apply (*it2, *it2e), ++ it2, ++ it2e));
00542 #endif
00543             ++ it1, ++ it1e;
00544         }
00545     }
00546     // Explicitly iterating column major
00547     template<template <class T1, class T2> class F, class M, class E>
00548     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00549     void iterating_matrix_assign (M &m, const matrix_expression<E> &e, column_major_tag) {
00550         typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
00551         typedef typename M::difference_type difference_type;
00552         difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
00553         difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
00554         typename M::iterator2 it2 (m.begin2 ());
00555         BOOST_UBLAS_CHECK (size1 == 0 || m.end2 () - it2 == size2, bad_size ());
00556         typename E::const_iterator2 it2e (e ().begin2 ());
00557         BOOST_UBLAS_CHECK (size1 == 0 || e ().end2 () - it2e == size2, bad_size ());
00558         while (-- size2 >= 0) {
00559 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00560             typename M::iterator1 it1 (it2.begin ());
00561             typename E::const_iterator1 it1e (it2e.begin ());
00562 #else
00563             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
00564             typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
00565 #endif
00566             BOOST_UBLAS_CHECK (it2.end () - it1 == size1, bad_size ());
00567             BOOST_UBLAS_CHECK (it2e.end () - it1e == size1, bad_size ());
00568             difference_type temp_size1 (size1);
00569 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
00570             while (-- temp_size1 >= 0)
00571                 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
00572 #else
00573             DD (temp_size1, 2, r, (functor_type::apply (*it1, *it1e), ++ it1, ++ it1e));
00574 #endif
00575             ++ it2, ++ it2e;
00576         }
00577     }
00578     // Explicitly indexing row major
00579     template<template <class T1, class T2> class F, class M, class E>
00580     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00581     void indexing_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
00582         typedef F<typename M::reference, typename E::value_type> functor_type;
00583         typedef typename M::size_type size_type;
00584         size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
00585         size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
00586         for (size_type i = 0; i < size1; ++ i) {
00587 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
00588             for (size_type j = 0; j < size2; ++ j)
00589                 functor_type::apply (m (i, j), e () (i, j));
00590 #else
00591             size_type j (0);
00592             DD (size2, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ j));
00593 #endif
00594         }
00595     }
00596     // Explicitly indexing column major
00597     template<template <class T1, class T2> class F, class M, class E>
00598     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00599     void indexing_matrix_assign (M &m, const matrix_expression<E> &e, column_major_tag) {
00600         typedef F<typename M::reference, typename E::value_type> functor_type;
00601         typedef typename M::size_type size_type;
00602         size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
00603         size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
00604         for (size_type j = 0; j < size2; ++ j) {
00605 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
00606             for (size_type i = 0; i < size1; ++ i)
00607                 functor_type::apply (m (i, j), e () (i, j));
00608 #else
00609             size_type i (0);
00610             DD (size1, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ i));
00611 #endif
00612         }
00613     }
00614 
00615     // Dense (proxy) case
00616     template<template <class T1, class T2> class F, class R, class M, class E, class C>
00617     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00618     void matrix_assign (M &m, const matrix_expression<E> &e, dense_proxy_tag, C) {
00619         // R unnecessary, make_conformant not required
00620         typedef C orientation_category;
00621 #ifdef BOOST_UBLAS_USE_INDEXING
00622         indexing_matrix_assign<F> (m, e, orientation_category ());
00623 #elif BOOST_UBLAS_USE_ITERATING
00624         iterating_matrix_assign<F> (m, e, orientation_category ());
00625 #else
00626         typedef typename M::difference_type difference_type;
00627         size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
00628         size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
00629         if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
00630             size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
00631             iterating_matrix_assign<F> (m, e, orientation_category ());
00632         else
00633             indexing_matrix_assign<F> (m, e, orientation_category ());
00634 #endif
00635     }
00636     // Packed (proxy) row major case
00637     template<template <class T1, class T2> class F, class R, class M, class E>
00638     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00639     void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
00640         typedef typename matrix_traits<E>::value_type expr_value_type;
00641         typedef F<typename M::iterator2::reference, expr_value_type> functor_type;
00642         // R unnecessary, make_conformant not required
00643         typedef typename M::difference_type difference_type;
00644         typedef typename M::value_type value_type;
00645         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
00646         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
00647 #if BOOST_UBLAS_TYPE_CHECK
00648         matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
00649         indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ());
00650         indexing_matrix_assign<F> (cm, e, row_major_tag ());
00651 #endif
00652         typename M::iterator1 it1 (m.begin1 ());
00653         typename M::iterator1 it1_end (m.end1 ());
00654         typename E::const_iterator1 it1e (e ().begin1 ());
00655         typename E::const_iterator1 it1e_end (e ().end1 ());
00656         difference_type it1_size (it1_end - it1);
00657         difference_type it1e_size (it1e_end - it1e);
00658         difference_type diff1 (0);
00659         if (it1_size > 0 && it1e_size > 0)
00660             diff1 = it1.index1 () - it1e.index1 ();
00661         if (diff1 != 0) {
00662             difference_type size1 = (std::min) (diff1, it1e_size);
00663             if (size1 > 0) {
00664                 it1e += size1;
00665                 it1e_size -= size1;
00666                 diff1 -= size1;
00667             }
00668             size1 = (std::min) (- diff1, it1_size);
00669             if (size1 > 0) {
00670                 it1_size -= size1;
00671                 if (!functor_type::computed) {
00672                     while (-- size1 >= 0) { // zeroing
00673 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00674                         typename M::iterator2 it2 (it1.begin ());
00675                         typename M::iterator2 it2_end (it1.end ());
00676 #else
00677                         typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00678                         typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
00679 #endif
00680                         difference_type size2 (it2_end - it2);
00681                         while (-- size2 >= 0)
00682                             functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
00683                         ++ it1;
00684                     }
00685                 } else {
00686                     it1 += size1;
00687                 }
00688                 diff1 += size1;
00689             }
00690         }
00691         difference_type size1 ((std::min) (it1_size, it1e_size));
00692         it1_size -= size1;
00693         it1e_size -= size1;
00694         while (-- size1 >= 0) {
00695 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00696             typename M::iterator2 it2 (it1.begin ());
00697             typename M::iterator2 it2_end (it1.end ());
00698             typename E::const_iterator2 it2e (it1e.begin ());
00699             typename E::const_iterator2 it2e_end (it1e.end ());
00700 #else
00701             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00702             typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
00703             typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
00704             typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
00705 #endif
00706             difference_type it2_size (it2_end - it2);
00707             difference_type it2e_size (it2e_end - it2e);
00708             difference_type diff2 (0);
00709             if (it2_size > 0 && it2e_size > 0) {
00710                 diff2 = it2.index2 () - it2e.index2 ();
00711                 difference_type size2 = (std::min) (diff2, it2e_size);
00712                 if (size2 > 0) {
00713                     it2e += size2;
00714                     it2e_size -= size2;
00715                     diff2 -= size2;
00716                 }
00717                 size2 = (std::min) (- diff2, it2_size);
00718                 if (size2 > 0) {
00719                     it2_size -= size2;
00720                     if (!functor_type::computed) {
00721                         while (-- size2 >= 0)   // zeroing
00722                             functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
00723                     } else {
00724                         it2 += size2;
00725                     }
00726                     diff2 += size2;
00727                 }
00728             }
00729             difference_type size2 ((std::min) (it2_size, it2e_size));
00730             it2_size -= size2;
00731             it2e_size -= size2;
00732             while (-- size2 >= 0)
00733                 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
00734             size2 = it2_size;
00735             if (!functor_type::computed) {
00736                 while (-- size2 >= 0)   // zeroing
00737                     functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
00738             } else {
00739                 it2 += size2;
00740             }
00741             ++ it1, ++ it1e;
00742         }
00743         size1 = it1_size;
00744         if (!functor_type::computed) {
00745             while (-- size1 >= 0) { // zeroing
00746 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00747                 typename M::iterator2 it2 (it1.begin ());
00748                 typename M::iterator2 it2_end (it1.end ());
00749 #else
00750                 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00751                 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
00752 #endif
00753                 difference_type size2 (it2_end - it2);
00754                 while (-- size2 >= 0)
00755                     functor_type::apply (*it2, expr_value_type/*zero*/()), ++ it2;
00756                 ++ it1;
00757             }
00758         } else {
00759             it1 += size1;
00760         }
00761 #if BOOST_UBLAS_TYPE_CHECK
00762         if (! disable_type_check<bool>::value)
00763             BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
00764 #endif
00765     }
00766     // Packed (proxy) column major case
00767     template<template <class T1, class T2> class F, class R, class M, class E>
00768     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00769     void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, column_major_tag) {
00770         typedef typename matrix_traits<E>::value_type expr_value_type;
00771         typedef F<typename M::iterator1::reference, expr_value_type> functor_type;
00772         // R unnecessary, make_conformant not required
00773         typedef typename M::difference_type difference_type;
00774         typedef typename M::value_type value_type;
00775         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
00776         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
00777 #if BOOST_UBLAS_TYPE_CHECK
00778         matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
00779         indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ());
00780         indexing_matrix_assign<F> (cm, e, column_major_tag ());
00781 #endif
00782         typename M::iterator2 it2 (m.begin2 ());
00783         typename M::iterator2 it2_end (m.end2 ());
00784         typename E::const_iterator2 it2e (e ().begin2 ());
00785         typename E::const_iterator2 it2e_end (e ().end2 ());
00786         difference_type it2_size (it2_end - it2);
00787         difference_type it2e_size (it2e_end - it2e);
00788         difference_type diff2 (0);
00789         if (it2_size > 0 && it2e_size > 0)
00790             diff2 = it2.index2 () - it2e.index2 ();
00791         if (diff2 != 0) {
00792             difference_type size2 = (std::min) (diff2, it2e_size);
00793             if (size2 > 0) {
00794                 it2e += size2;
00795                 it2e_size -= size2;
00796                 diff2 -= size2;
00797             }
00798             size2 = (std::min) (- diff2, it2_size);
00799             if (size2 > 0) {
00800                 it2_size -= size2;
00801                 if (!functor_type::computed) {
00802                     while (-- size2 >= 0) { // zeroing
00803 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00804                         typename M::iterator1 it1 (it2.begin ());
00805                         typename M::iterator1 it1_end (it2.end ());
00806 #else
00807                         typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
00808                         typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
00809 #endif
00810                         difference_type size1 (it1_end - it1);
00811                         while (-- size1 >= 0)
00812                             functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
00813                         ++ it2;
00814                     }
00815                 } else {
00816                     it2 += size2;
00817                 }
00818                 diff2 += size2;
00819             }
00820         }
00821         difference_type size2 ((std::min) (it2_size, it2e_size));
00822         it2_size -= size2;
00823         it2e_size -= size2;
00824         while (-- size2 >= 0) {
00825 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00826             typename M::iterator1 it1 (it2.begin ());
00827             typename M::iterator1 it1_end (it2.end ());
00828             typename E::const_iterator1 it1e (it2e.begin ());
00829             typename E::const_iterator1 it1e_end (it2e.end ());
00830 #else
00831             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
00832             typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
00833             typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
00834             typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
00835 #endif
00836             difference_type it1_size (it1_end - it1);
00837             difference_type it1e_size (it1e_end - it1e);
00838             difference_type diff1 (0);
00839             if (it1_size > 0 && it1e_size > 0) {
00840                 diff1 = it1.index1 () - it1e.index1 ();
00841                 difference_type size1 = (std::min) (diff1, it1e_size);
00842                 if (size1 > 0) {
00843                     it1e += size1;
00844                     it1e_size -= size1;
00845                     diff1 -= size1;
00846                 }
00847                 size1 = (std::min) (- diff1, it1_size);
00848                 if (size1 > 0) {
00849                     it1_size -= size1;
00850                     if (!functor_type::computed) {
00851                         while (-- size1 >= 0)   // zeroing
00852                             functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
00853                     } else {
00854                         it1 += size1;
00855                     }
00856                     diff1 += size1;
00857                 }
00858             }
00859             difference_type size1 ((std::min) (it1_size, it1e_size));
00860             it1_size -= size1;
00861             it1e_size -= size1;
00862             while (-- size1 >= 0)
00863                 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
00864             size1 = it1_size;
00865             if (!functor_type::computed) {
00866                 while (-- size1 >= 0)   // zeroing
00867                     functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
00868             } else {
00869                 it1 += size1;
00870             }
00871             ++ it2, ++ it2e;
00872         }
00873         size2 = it2_size;
00874         if (!functor_type::computed) {
00875             while (-- size2 >= 0) { // zeroing
00876 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00877                 typename M::iterator1 it1 (it2.begin ());
00878                 typename M::iterator1 it1_end (it2.end ());
00879 #else
00880                 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
00881                 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
00882 #endif
00883                 difference_type size1 (it1_end - it1);
00884                 while (-- size1 >= 0)
00885                     functor_type::apply (*it1, expr_value_type/*zero*/()), ++ it1;
00886                 ++ it2;
00887             }
00888         } else {
00889             it2 += size2;
00890         }
00891 #if BOOST_UBLAS_TYPE_CHECK
00892         if (! disable_type_check<bool>::value)
00893             BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
00894 #endif
00895     }
00896     // Sparse row major case
00897     template<template <class T1, class T2> class F, class R, class M, class E>
00898     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00899     void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, row_major_tag) {
00900         typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
00901         // R unnecessary, make_conformant not required
00902         BOOST_STATIC_ASSERT ((!functor_type::computed));
00903         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
00904         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
00905         typedef typename M::value_type value_type;
00906         // Sparse type has no numeric constraints to check
00907 
00908         m.clear ();
00909         typename E::const_iterator1 it1e (e ().begin1 ());
00910         typename E::const_iterator1 it1e_end (e ().end1 ());
00911         while (it1e != it1e_end) {
00912 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00913             typename E::const_iterator2 it2e (it1e.begin ());
00914             typename E::const_iterator2 it2e_end (it1e.end ());
00915 #else
00916             typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
00917             typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
00918 #endif
00919             while (it2e != it2e_end) {
00920                 value_type t (*it2e);
00921                 if (t != value_type/*zero*/())
00922                     m.insert_element (it2e.index1 (), it2e.index2 (), t);
00923                 ++ it2e;
00924             }
00925             ++ it1e;
00926         }
00927     }
00928     // Sparse column major case
00929     template<template <class T1, class T2> class F, class R, class M, class E>
00930     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00931     void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, column_major_tag) {
00932         typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
00933         // R unnecessary, make_conformant not required
00934         BOOST_STATIC_ASSERT ((!functor_type::computed));
00935         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
00936         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
00937         typedef typename M::value_type value_type;
00938         // Sparse type has no numeric constraints to check
00939 
00940         m.clear ();
00941         typename E::const_iterator2 it2e (e ().begin2 ());
00942         typename E::const_iterator2 it2e_end (e ().end2 ());
00943         while (it2e != it2e_end) {
00944 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00945             typename E::const_iterator1 it1e (it2e.begin ());
00946             typename E::const_iterator1 it1e_end (it2e.end ());
00947 #else
00948             typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
00949             typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
00950 #endif
00951             while (it1e != it1e_end) {
00952                 value_type t (*it1e);
00953                 if (t != value_type/*zero*/())
00954                     m.insert_element (it1e.index1 (), it1e.index2 (), t);
00955                 ++ it1e;
00956             }
00957             ++ it2e;
00958         }
00959     }
00960     // Sparse proxy or functional row major case
00961     template<template <class T1, class T2> class F, class R, class M, class E>
00962     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
00963     void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) {
00964         typedef typename matrix_traits<E>::value_type expr_value_type;
00965         typedef F<typename M::iterator2::reference, expr_value_type> functor_type;
00966         typedef R conformant_restrict_type;
00967         typedef typename M::size_type size_type;
00968         typedef typename M::difference_type difference_type;
00969         typedef typename M::value_type value_type;
00970         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
00971         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
00972 #if BOOST_UBLAS_TYPE_CHECK
00973         matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
00974         indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ());
00975         indexing_matrix_assign<F> (cm, e, row_major_tag ());
00976 #endif
00977         detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
00978 
00979         typename M::iterator1 it1 (m.begin1 ());
00980         typename M::iterator1 it1_end (m.end1 ());
00981         typename E::const_iterator1 it1e (e ().begin1 ());
00982         typename E::const_iterator1 it1e_end (e ().end1 ());
00983         while (it1 != it1_end && it1e != it1e_end) {
00984             difference_type compare = it1.index1 () - it1e.index1 ();
00985             if (compare == 0) {
00986 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00987                 typename M::iterator2 it2 (it1.begin ());
00988                 typename M::iterator2 it2_end (it1.end ());
00989                 typename E::const_iterator2 it2e (it1e.begin ());
00990                 typename E::const_iterator2 it2e_end (it1e.end ());
00991 #else
00992                 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
00993                 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
00994                 typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
00995                 typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
00996 #endif
00997                 if (it2 != it2_end && it2e != it2e_end) {
00998                     size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
00999                     while (true) {
01000                         difference_type compare = it2_index - it2e_index;
01001                         if (compare == 0) {
01002                             functor_type::apply (*it2, *it2e);
01003                             ++ it2, ++ it2e;
01004                             if (it2 != it2_end && it2e != it2e_end) {
01005                                 it2_index = it2.index2 ();
01006                                 it2e_index = it2e.index2 ();
01007                             } else
01008                                 break;
01009                         } else if (compare < 0) {
01010                             if (!functor_type::computed) {
01011                                 functor_type::apply (*it2, expr_value_type/*zero*/());
01012                                 ++ it2;
01013                             } else
01014                                 increment (it2, it2_end, - compare);
01015                             if (it2 != it2_end)
01016                                 it2_index = it2.index2 ();
01017                             else
01018                                 break;
01019                         } else if (compare > 0) {
01020                             increment (it2e, it2e_end, compare);
01021                             if (it2e != it2e_end)
01022                                 it2e_index = it2e.index2 ();
01023                             else
01024                                 break;
01025                         }
01026                     }
01027                 }
01028                 if (!functor_type::computed) {
01029                     while (it2 != it2_end) {    // zeroing
01030                         functor_type::apply (*it2, expr_value_type/*zero*/());
01031                         ++ it2;
01032                     }
01033                 } else {
01034                     it2 = it2_end;
01035                 }
01036                 ++ it1, ++ it1e;
01037             } else if (compare < 0) {
01038                 if (!functor_type::computed) {
01039 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01040                     typename M::iterator2 it2 (it1.begin ());
01041                     typename M::iterator2 it2_end (it1.end ());
01042 #else
01043                     typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
01044                     typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
01045 #endif
01046                     while (it2 != it2_end) {    // zeroing
01047                         functor_type::apply (*it2, expr_value_type/*zero*/());
01048                         ++ it2;
01049                     }
01050                     ++ it1;
01051                 } else {
01052                     increment (it1, it1_end, - compare);
01053                 }
01054             } else if (compare > 0) {
01055                 increment (it1e, it1e_end, compare);
01056             }
01057         }
01058         if (!functor_type::computed) {
01059             while (it1 != it1_end) {
01060 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01061                 typename M::iterator2 it2 (it1.begin ());
01062                 typename M::iterator2 it2_end (it1.end ());
01063 #else
01064                 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
01065                 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
01066 #endif
01067                 while (it2 != it2_end) {    // zeroing
01068                     functor_type::apply (*it2, expr_value_type/*zero*/());
01069                     ++ it2;
01070                 }
01071                 ++ it1;
01072             }
01073         } else {
01074             it1 = it1_end;
01075         }
01076 #if BOOST_UBLAS_TYPE_CHECK
01077         if (! disable_type_check<bool>::value)
01078             BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
01079 #endif
01080     }
01081     // Sparse proxy or functional column major case
01082     template<template <class T1, class T2> class F, class R, class M, class E>
01083     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
01084     void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) {
01085         typedef typename matrix_traits<E>::value_type expr_value_type;
01086         typedef F<typename M::iterator1::reference, expr_value_type> functor_type;
01087         typedef R conformant_restrict_type;
01088         typedef typename M::size_type size_type;
01089         typedef typename M::difference_type difference_type;
01090         typedef typename M::value_type value_type;
01091         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
01092         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
01093 #if BOOST_UBLAS_TYPE_CHECK
01094         matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
01095         indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ());
01096         indexing_matrix_assign<F> (cm, e, column_major_tag ());
01097 #endif
01098         detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
01099 
01100         typename M::iterator2 it2 (m.begin2 ());
01101         typename M::iterator2 it2_end (m.end2 ());
01102         typename E::const_iterator2 it2e (e ().begin2 ());
01103         typename E::const_iterator2 it2e_end (e ().end2 ());
01104         while (it2 != it2_end && it2e != it2e_end) {
01105             difference_type compare = it2.index2 () - it2e.index2 ();
01106             if (compare == 0) {
01107 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01108                 typename M::iterator1 it1 (it2.begin ());
01109                 typename M::iterator1 it1_end (it2.end ());
01110                 typename E::const_iterator1 it1e (it2e.begin ());
01111                 typename E::const_iterator1 it1e_end (it2e.end ());
01112 #else
01113                 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
01114                 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
01115                 typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
01116                 typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
01117 #endif
01118                 if (it1 != it1_end && it1e != it1e_end) {
01119                     size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
01120                     while (true) {
01121                         difference_type compare = it1_index - it1e_index;
01122                         if (compare == 0) {
01123                             functor_type::apply (*it1, *it1e);
01124                             ++ it1, ++ it1e;
01125                             if (it1 != it1_end && it1e != it1e_end) {
01126                                 it1_index = it1.index1 ();
01127                                 it1e_index = it1e.index1 ();
01128                             } else
01129                                 break;
01130                         } else if (compare < 0) {
01131                             if (!functor_type::computed) {
01132                                 functor_type::apply (*it1, expr_value_type/*zero*/()); // zeroing
01133                                 ++ it1;
01134                             } else
01135                                 increment (it1, it1_end, - compare);
01136                             if (it1 != it1_end)
01137                                 it1_index = it1.index1 ();
01138                             else
01139                                 break;
01140                         } else if (compare > 0) {
01141                             increment (it1e, it1e_end, compare);
01142                             if (it1e != it1e_end)
01143                                 it1e_index = it1e.index1 ();
01144                             else
01145                                 break;
01146                         }
01147                     }
01148                 }
01149                 if (!functor_type::computed) {
01150                     while (it1 != it1_end) {    // zeroing
01151                         functor_type::apply (*it1, expr_value_type/*zero*/());
01152                         ++ it1;
01153                     }
01154                 } else {
01155                     it1 = it1_end;
01156                 }
01157                 ++ it2, ++ it2e;
01158             } else if (compare < 0) {
01159                 if (!functor_type::computed) {
01160 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01161                     typename M::iterator1 it1 (it2.begin ());
01162                     typename M::iterator1 it1_end (it2.end ());
01163 #else
01164                     typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
01165                     typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
01166 #endif
01167                     while (it1 != it1_end) {    // zeroing
01168                         functor_type::apply (*it1, expr_value_type/*zero*/());
01169                         ++ it1;
01170                     }
01171                     ++ it2;
01172                 } else {
01173                     increment (it2, it2_end, - compare);
01174                 }
01175             } else if (compare > 0) {
01176                 increment (it2e, it2e_end, compare);
01177             }
01178         }
01179         if (!functor_type::computed) {
01180             while (it2 != it2_end) {
01181 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01182                 typename M::iterator1 it1 (it2.begin ());
01183                 typename M::iterator1 it1_end (it2.end ());
01184 #else
01185                 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
01186                 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
01187 #endif
01188                 while (it1 != it1_end) {    // zeroing
01189                     functor_type::apply (*it1, expr_value_type/*zero*/());
01190                     ++ it1;
01191                 }
01192                 ++ it2;
01193             }
01194         } else {
01195             it2 = it2_end;
01196         }
01197 #if BOOST_UBLAS_TYPE_CHECK
01198         if (! disable_type_check<bool>::value)
01199             BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
01200 #endif
01201     }
01202 
01203     // Dispatcher
01204     template<template <class T1, class T2> class F, class M, class E>
01205     BOOST_UBLAS_INLINE
01206     void matrix_assign (M &m, const matrix_expression<E> &e) {
01207         typedef typename matrix_assign_traits<typename M::storage_category,
01208                                               F<typename M::reference, typename E::value_type>::computed,
01209                                               typename E::const_iterator1::iterator_category,
01210                                               typename E::const_iterator2::iterator_category>::storage_category storage_category;
01211         // give preference to matrix M's orientation if known
01212         typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
01213                                           typename E::orientation_category ,
01214                                           typename M::orientation_category >::type orientation_category;
01215         typedef basic_full<typename M::size_type> unrestricted;
01216         matrix_assign<F, unrestricted> (m, e, storage_category (), orientation_category ());
01217     }
01218     template<template <class T1, class T2> class F, class R, class M, class E>
01219     BOOST_UBLAS_INLINE
01220     void matrix_assign (M &m, const matrix_expression<E> &e) {
01221         typedef R conformant_restrict_type;
01222         typedef typename matrix_assign_traits<typename M::storage_category,
01223                                               F<typename M::reference, typename E::value_type>::computed,
01224                                               typename E::const_iterator1::iterator_category,
01225                                               typename E::const_iterator2::iterator_category>::storage_category storage_category;
01226         // give preference to matrix M's orientation if known
01227         typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
01228                                           typename E::orientation_category ,
01229                                           typename M::orientation_category >::type orientation_category;
01230         matrix_assign<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ());
01231     }
01232 
01233     template<class SC, class RI1, class RI2>
01234     struct matrix_swap_traits {
01235         typedef SC storage_category;
01236     };
01237 
01238     template<>
01239     struct matrix_swap_traits<dense_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
01240         typedef sparse_proxy_tag storage_category;
01241     };
01242 
01243     template<>
01244     struct matrix_swap_traits<packed_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
01245         typedef sparse_proxy_tag storage_category;
01246     };
01247 
01248     // Dense (proxy) row major case
01249     template<template <class T1, class T2> class F, class R, class M, class E>
01250     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
01251     void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, row_major_tag) {
01252         typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
01253         // R unnecessary, make_conformant not required
01254         typedef typename M::size_type size_type;
01255         typedef typename M::difference_type difference_type;
01256         typename M::iterator1 it1 (m.begin1 ());
01257         typename E::iterator1 it1e (e ().begin1 ());
01258         difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (e ().end1 () - it1e)));
01259         while (-- size1 >= 0) {
01260 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01261             typename M::iterator2 it2 (it1.begin ());
01262             typename E::iterator2 it2e (it1e.begin ());
01263             difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (it1e.end () - it2e)));
01264 #else
01265             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
01266             typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
01267             difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (end (it1e, iterator1_tag ()) - it2e)));
01268 #endif
01269             while (-- size2 >= 0)
01270                 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
01271             ++ it1, ++ it1e;
01272         }
01273     }
01274     // Dense (proxy) column major case
01275     template<template <class T1, class T2> class F, class R, class M, class E>
01276     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
01277     void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, column_major_tag) {
01278         typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
01279         // R unnecessary, make_conformant not required
01280         typedef typename M::size_type size_type;
01281         typedef typename M::difference_type difference_type;
01282         typename M::iterator2 it2 (m.begin2 ());
01283         typename E::iterator2 it2e (e ().begin2 ());
01284         difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (e ().end2 () - it2e)));
01285         while (-- size2 >= 0) {
01286 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01287             typename M::iterator1 it1 (it2.begin ());
01288             typename E::iterator1 it1e (it2e.begin ());
01289             difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (it2e.end () - it1e)));
01290 #else
01291             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
01292             typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
01293             difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (end (it2e, iterator2_tag ()) - it1e)));
01294 #endif
01295             while (-- size1 >= 0)
01296                 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
01297             ++ it2, ++ it2e;
01298         }
01299     }
01300     // Packed (proxy) row major case
01301     template<template <class T1, class T2> class F, class R, class M, class E>
01302     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
01303     void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
01304         typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
01305         // R unnecessary, make_conformant not required
01306         typedef typename M::size_type size_type;
01307         typedef typename M::difference_type difference_type;
01308         typename M::iterator1 it1 (m.begin1 ());
01309         typename E::iterator1 it1e (e ().begin1 ());
01310         difference_type size1 (BOOST_UBLAS_SAME (m.end1 () - it1, e ().end1 () - it1e));
01311         while (-- size1 >= 0) {
01312 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01313             typename M::iterator2 it2 (it1.begin ());
01314             typename E::iterator2 it2e (it1e.begin ());
01315             difference_type size2 (BOOST_UBLAS_SAME (it1.end () - it2, it1e.end () - it2e));
01316 #else
01317             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
01318             typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
01319             difference_type size2 (BOOST_UBLAS_SAME (end (it1, iterator1_tag ()) - it2, end (it1e, iterator1_tag ()) - it2e));
01320 #endif
01321             while (-- size2 >= 0)
01322                 functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
01323             ++ it1, ++ it1e;
01324         }
01325     }
01326     // Packed (proxy) column major case
01327     template<template <class T1, class T2> class F, class R, class M, class E>
01328     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
01329     void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, column_major_tag) {
01330         typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
01331         // R unnecessary, make_conformant not required
01332         typedef typename M::size_type size_type;
01333         typedef typename M::difference_type difference_type;
01334         typename M::iterator2 it2 (m.begin2 ());
01335         typename E::iterator2 it2e (e ().begin2 ());
01336         difference_type size2 (BOOST_UBLAS_SAME (m.end2 () - it2, e ().end2 () - it2e));
01337         while (-- size2 >= 0) {
01338 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01339             typename M::iterator1 it1 (it2.begin ());
01340             typename E::iterator1 it1e (it2e.begin ());
01341             difference_type size1 (BOOST_UBLAS_SAME (it2.end () - it1, it2e.end () - it1e));
01342 #else
01343             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
01344             typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
01345             difference_type size1 (BOOST_UBLAS_SAME (end (it2, iterator2_tag ()) - it1, end (it2e, iterator2_tag ()) - it1e));
01346 #endif
01347             while (-- size1 >= 0)
01348                 functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
01349             ++ it2, ++ it2e;
01350         }
01351     }
01352     // Sparse (proxy) row major case
01353     template<template <class T1, class T2> class F, class R, class M, class E>
01354     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
01355     void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) {
01356         typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
01357         typedef R conformant_restrict_type;
01358         typedef typename M::size_type size_type;
01359         typedef typename M::difference_type difference_type;
01360         typedef typename M::value_type value_type;
01361         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
01362         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
01363 
01364         detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
01365         // FIXME should be a seperate restriction for E
01366         detail::make_conformant (e (), m, row_major_tag (), conformant_restrict_type ());
01367 
01368         typename M::iterator1 it1 (m.begin1 ());
01369         typename M::iterator1 it1_end (m.end1 ());
01370         typename E::iterator1 it1e (e ().begin1 ());
01371         typename E::iterator1 it1e_end (e ().end1 ());
01372         while (it1 != it1_end && it1e != it1e_end) {
01373             difference_type compare = it1.index1 () - it1e.index1 ();
01374             if (compare == 0) {
01375 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01376                 typename M::iterator2 it2 (it1.begin ());
01377                 typename M::iterator2 it2_end (it1.end ());
01378                 typename E::iterator2 it2e (it1e.begin ());
01379                 typename E::iterator2 it2e_end (it1e.end ());
01380 #else
01381                 typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
01382                 typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
01383                 typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
01384                 typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
01385 #endif
01386                 if (it2 != it2_end && it2e != it2e_end) {
01387                     size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
01388                     while (true) {
01389                         difference_type compare = it2_index - it2e_index;
01390                         if (compare == 0) {
01391                             functor_type::apply (*it2, *it2e);
01392                             ++ it2, ++ it2e;
01393                             if (it2 != it2_end && it2e != it2e_end) {
01394                                 it2_index = it2.index2 ();
01395                                 it2e_index = it2e.index2 ();
01396                             } else
01397                                 break;
01398                         } else if (compare < 0) {
01399                             increment (it2, it2_end, - compare);
01400                             if (it2 != it2_end)
01401                                 it2_index = it2.index2 ();
01402                             else
01403                                 break;
01404                         } else if (compare > 0) {
01405                             increment (it2e, it2e_end, compare);
01406                             if (it2e != it2e_end)
01407                                 it2e_index = it2e.index2 ();
01408                             else
01409                                 break;
01410                         }
01411                     }
01412                 }
01413 #if BOOST_UBLAS_TYPE_CHECK
01414                 increment (it2e, it2e_end);
01415                 increment (it2, it2_end);
01416 #endif
01417                 ++ it1, ++ it1e;
01418             } else if (compare < 0) {
01419 #if BOOST_UBLAS_TYPE_CHECK
01420                 while (it1.index1 () < it1e.index1 ()) {
01421 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01422                     typename M::iterator2 it2 (it1.begin ());
01423                     typename M::iterator2 it2_end (it1.end ());
01424 #else
01425                     typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
01426                     typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
01427 #endif
01428                     increment (it2, it2_end);
01429                     ++ it1;
01430                 }
01431 #else
01432                 increment (it1, it1_end, - compare);
01433 #endif
01434             } else if (compare > 0) {
01435 #if BOOST_UBLAS_TYPE_CHECK
01436                 while (it1e.index1 () < it1.index1 ()) {
01437 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01438                     typename E::iterator2 it2e (it1e.begin ());
01439                     typename E::iterator2 it2e_end (it1e.end ());
01440 #else
01441                     typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
01442                     typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
01443 #endif
01444                     increment (it2e, it2e_end);
01445                     ++ it1e;
01446                 }
01447 #else
01448                 increment (it1e, it1e_end, compare);
01449 #endif
01450             }
01451         }
01452 #if BOOST_UBLAS_TYPE_CHECK
01453         while (it1e != it1e_end) {
01454 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01455             typename E::iterator2 it2e (it1e.begin ());
01456             typename E::iterator2 it2e_end (it1e.end ());
01457 #else
01458             typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
01459             typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
01460 #endif
01461             increment (it2e, it2e_end);
01462             ++ it1e;
01463         }
01464         while (it1 != it1_end) {
01465 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01466             typename M::iterator2 it2 (it1.begin ());
01467             typename M::iterator2 it2_end (it1.end ());
01468 #else
01469             typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
01470             typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
01471 #endif
01472             increment (it2, it2_end);
01473             ++ it1;
01474         }
01475 #endif
01476     }
01477     // Sparse (proxy) column major case
01478     template<template <class T1, class T2> class F, class R, class M, class E>
01479     // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
01480     void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) {
01481         typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
01482         typedef R conformant_restrict_type;
01483         typedef typename M::size_type size_type;
01484         typedef typename M::difference_type difference_type;
01485         typedef typename M::value_type value_type;
01486         BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
01487         BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
01488 
01489         detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
01490         // FIXME should be a seperate restriction for E
01491         detail::make_conformant (e (), m, column_major_tag (), conformant_restrict_type ());
01492 
01493         typename M::iterator2 it2 (m.begin2 ());
01494         typename M::iterator2 it2_end (m.end2 ());
01495         typename E::iterator2 it2e (e ().begin2 ());
01496         typename E::iterator2 it2e_end (e ().end2 ());
01497         while (it2 != it2_end && it2e != it2e_end) {
01498             difference_type compare = it2.index2 () - it2e.index2 ();
01499             if (compare == 0) {
01500 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01501                 typename M::iterator1 it1 (it2.begin ());
01502                 typename M::iterator1 it1_end (it2.end ());
01503                 typename E::iterator1 it1e (it2e.begin ());
01504                 typename E::iterator1 it1e_end (it2e.end ());
01505 #else
01506                 typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
01507                 typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
01508                 typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
01509                 typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
01510 #endif
01511                 if (it1 != it1_end && it1e != it1e_end) {
01512                     size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
01513                     while (true) {
01514                         difference_type compare = it1_index - it1e_index;
01515                         if (compare == 0) {
01516                             functor_type::apply (*it1, *it1e);
01517                             ++ it1, ++ it1e;
01518                             if (it1 != it1_end && it1e != it1e_end) {
01519                                 it1_index = it1.index1 ();
01520                                 it1e_index = it1e.index1 ();
01521                             } else
01522                                 break;
01523                         }  else if (compare < 0) {
01524                             increment (it1, it1_end, - compare);
01525                             if (it1 != it1_end)
01526                                 it1_index = it1.index1 ();
01527                             else
01528                                 break;
01529                         } else if (compare > 0) {
01530                             increment (it1e, it1e_end, compare);
01531                             if (it1e != it1e_end)
01532                                 it1e_index = it1e.index1 ();
01533                             else
01534                                 break;
01535                         }
01536                     }
01537                 }
01538 #if BOOST_UBLAS_TYPE_CHECK
01539                 increment (it1e, it1e_end);
01540                 increment (it1, it1_end);
01541 #endif
01542                 ++ it2, ++ it2e;
01543             } else if (compare < 0) {
01544 #if BOOST_UBLAS_TYPE_CHECK
01545                 while (it2.index2 () < it2e.index2 ()) {
01546 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01547                     typename M::iterator1 it1 (it2.begin ());
01548                     typename M::iterator1 it1_end (it2.end ());
01549 #else
01550                     typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
01551                     typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
01552 #endif
01553                     increment (it1, it1_end);
01554                     ++ it2;
01555                 }
01556 #else
01557                 increment (it2, it2_end, - compare);
01558 #endif
01559             } else if (compare > 0) {
01560 #if BOOST_UBLAS_TYPE_CHECK
01561                 while (it2e.index2 () < it2.index2 ()) {
01562 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01563                     typename E::iterator1 it1e (it2e.begin ());
01564                     typename E::iterator1 it1e_end (it2e.end ());
01565 #else
01566                     typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
01567                     typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
01568 #endif
01569                     increment (it1e, it1e_end);
01570                     ++ it2e;
01571                 }
01572 #else
01573                 increment (it2e, it2e_end, compare);
01574 #endif
01575             }
01576         }
01577 #if BOOST_UBLAS_TYPE_CHECK
01578         while (it2e != it2e_end) {
01579 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01580             typename E::iterator1 it1e (it2e.begin ());
01581             typename E::iterator1 it1e_end (it2e.end ());
01582 #else
01583             typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
01584             typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
01585 #endif
01586             increment (it1e, it1e_end);
01587             ++ it2e;
01588         }
01589         while (it2 != it2_end) {
01590 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01591             typename M::iterator1 it1 (it2.begin ());
01592             typename M::iterator1 it1_end (it2.end ());
01593 #else
01594             typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
01595             typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
01596 #endif
01597             increment (it1, it1_end);
01598             ++ it2;
01599         }
01600 #endif
01601     }
01602 
01603     // Dispatcher
01604     template<template <class T1, class T2> class F, class M, class E>
01605     BOOST_UBLAS_INLINE
01606     void matrix_swap (M &m, matrix_expression<E> &e) {
01607         typedef typename matrix_swap_traits<typename M::storage_category,
01608                                             typename E::const_iterator1::iterator_category,
01609                                             typename E::const_iterator2::iterator_category>::storage_category storage_category;
01610         // give preference to matrix M's orientation if known
01611         typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
01612                                           typename E::orientation_category ,
01613                                           typename M::orientation_category >::type orientation_category;
01614         typedef basic_full<typename M::size_type> unrestricted;
01615         matrix_swap<F, unrestricted> (m, e, storage_category (), orientation_category ());
01616     }
01617     template<template <class T1, class T2> class F, class R, class M, class E>
01618     BOOST_UBLAS_INLINE
01619     void matrix_swap (M &m, matrix_expression<E> &e) {
01620         typedef R conformant_restrict_type;
01621         typedef typename matrix_swap_traits<typename M::storage_category,
01622                                             typename E::const_iterator1::iterator_category,
01623                                             typename E::const_iterator2::iterator_category>::storage_category storage_category;
01624         // give preference to matrix M's orientation if known
01625         typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
01626                                           typename E::orientation_category ,
01627                                           typename M::orientation_category >::type orientation_category;
01628         matrix_swap<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ());
01629     }
01630 
01631 }}}
01632 
01633 #endif