source: NonGTP/Boost/boost/numeric/ublas/detail/matrix_assign.hpp @ 857

Revision 857, 78.1 KB checked in by igarcia, 19 years ago (diff)
Line 
1//
2//  Copyright (c) 2000-2002
3//  Joerg Walter, Mathias Koch
4//
5//  Permission to use, copy, modify, distribute and sell this software
6//  and its documentation for any purpose is hereby granted without fee,
7//  provided that the above copyright notice appear in all copies and
8//  that both that copyright notice and this permission notice appear
9//  in supporting documentation.  The authors make no representations
10//  about the suitability of this software for any purpose.
11//  It is provided "as is" without express or implied warranty.
12//
13//  The authors gratefully acknowledge the support of
14//  GeNeSys mbH & Co. KG in producing this work.
15//
16
17#ifndef _BOOST_UBLAS_MATRIX_ASSIGN_
18#define _BOOST_UBLAS_MATRIX_ASSIGN_
19
20// Required for make_conformant storage
21#include <vector>
22
23// Iterators based on ideas of Jeremy Siek
24
25namespace boost { namespace numeric { namespace ublas {
26namespace detail {
27   
28    // Weak equality check - useful to compare equality two arbitary matrix expression results.
29    // Since the actual expressions are unknown, we check for and arbitary error bound
30    // on the relative error.
31    // For a linear expression the infinity norm makes sense as we do not know how the elements will be
32    // combined in the expression. False positive results are inevitable for arbirary expressions!
33    template<class E1, class E2, class S>
34    BOOST_UBLAS_INLINE
35    bool equals (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2, S epsilon, S min_norm) {
36        return norm_inf (e1 - e2) < epsilon *
37               std::max<S> (std::max<S> (norm_inf (e1), norm_inf (e2)), min_norm);
38    }
39
40    template<class E1, class E2>
41    BOOST_UBLAS_INLINE
42    bool expression_type_check (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2) {
43        typedef typename type_traits<typename promote_traits<typename E1::value_type,
44                                     typename E2::value_type>::promote_type>::real_type real_type;
45        return equals (e1, e2, BOOST_UBLAS_TYPE_CHECK_EPSILON, BOOST_UBLAS_TYPE_CHECK_MIN);
46    }
47
48
49    template<class M, class E, class R>
50    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
51    void make_conformant (M &m, const matrix_expression<E> &e, row_major_tag, R) {
52        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
53        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
54        typedef R conformant_restrict_type;
55        typedef typename M::size_type size_type;
56        typedef typename M::difference_type difference_type;
57        typedef typename M::value_type value_type;
58        // FIXME unbounded_array with push_back maybe better
59        std::vector<std::pair<size_type, size_type> > index;
60        typename M::iterator1 it1 (m.begin1 ());
61        typename M::iterator1 it1_end (m.end1 ());
62        typename E::const_iterator1 it1e (e ().begin1 ());
63        typename E::const_iterator1 it1e_end (e ().end1 ());
64        while (it1 != it1_end && it1e != it1e_end) {
65            difference_type compare = it1.index1 () - it1e.index1 ();
66            if (compare == 0) {
67#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
68                typename M::iterator2 it2 (it1.begin ());
69                typename M::iterator2 it2_end (it1.end ());
70                typename E::const_iterator2 it2e (it1e.begin ());
71                typename E::const_iterator2 it2e_end (it1e.end ());
72#else
73                typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
74                typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
75                typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
76                typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
77#endif
78                if (it2 != it2_end && it2e != it2e_end) {
79                    size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
80                    while (true) {
81                        difference_type compare = it2_index - it2e_index;
82                        if (compare == 0) {
83                            ++ it2, ++ it2e;
84                            if (it2 != it2_end && it2e != it2e_end) {
85                                it2_index = it2.index2 ();
86                                it2e_index = it2e.index2 ();
87                            } else
88                                break;
89                        } else if (compare < 0) {
90                            increment (it2, it2_end, - compare);
91                            if (it2 != it2_end)
92                                it2_index = it2.index2 ();
93                            else
94                                break;
95                        } else if (compare > 0) {
96                            if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
97                                if (*it2e != value_type/*zero*/())
98                                    index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
99                            ++ it2e;
100                            if (it2e != it2e_end)
101                                it2e_index = it2e.index2 ();
102                            else
103                                break;
104                        }
105                    }
106                }
107                while (it2e != it2e_end) {
108                    if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
109                        if (*it2e != value_type/*zero*/())
110                            index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
111                    ++ it2e;
112                }
113                ++ it1, ++ it1e;
114            } else if (compare < 0) {
115                increment (it1, it1_end, - compare);
116            } else if (compare > 0) {
117#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
118                typename E::const_iterator2 it2e (it1e.begin ());
119                typename E::const_iterator2 it2e_end (it1e.end ());
120#else
121                typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
122                typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
123#endif
124                while (it2e != it2e_end) {
125                    if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
126                        if (*it2e != value_type/*zero*/())
127                            index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
128                    ++ it2e;
129                }
130                ++ it1e;
131            }
132        }
133        while (it1e != it1e_end) {
134#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
135            typename E::const_iterator2 it2e (it1e.begin ());
136            typename E::const_iterator2 it2e_end (it1e.end ());
137#else
138            typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
139            typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
140#endif
141            while (it2e != it2e_end) {
142                if (conformant_restrict_type::other (it2e.index1 (), it2e.index2 ()))
143                    if (*it2e != value_type/*zero*/())
144                        index.push_back (std::pair<size_type, size_type> (it2e.index1 (), it2e.index2 ()));
145                ++ it2e;
146            }
147            ++ it1e;
148        }
149        // ISSUE proxies require insert_element
150        for (size_type k = 0; k < index.size (); ++ k)
151            m (index [k].first, index [k].second) = value_type/*zero*/();
152    }
153    template<class M, class E, class R>
154    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
155    void make_conformant (M &m, const matrix_expression<E> &e, column_major_tag, R) {
156        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
157        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
158        typedef R conformant_restrict_type;
159        typedef typename M::size_type size_type;
160        typedef typename M::difference_type difference_type;
161        typedef typename M::value_type value_type;
162        std::vector<std::pair<size_type, size_type> > index;
163        typename M::iterator2 it2 (m.begin2 ());
164        typename M::iterator2 it2_end (m.end2 ());
165        typename E::const_iterator2 it2e (e ().begin2 ());
166        typename E::const_iterator2 it2e_end (e ().end2 ());
167        while (it2 != it2_end && it2e != it2e_end) {
168            difference_type compare = it2.index2 () - it2e.index2 ();
169            if (compare == 0) {
170#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
171                typename M::iterator1 it1 (it2.begin ());
172                typename M::iterator1 it1_end (it2.end ());
173                typename E::const_iterator1 it1e (it2e.begin ());
174                typename E::const_iterator1 it1e_end (it2e.end ());
175#else
176                typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
177                typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
178                typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
179                typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
180#endif
181                if (it1 != it1_end && it1e != it1e_end) {
182                    size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
183                    while (true) {
184                        difference_type compare = it1_index - it1e_index;
185                        if (compare == 0) {
186                            ++ it1, ++ it1e;
187                            if (it1 != it1_end && it1e != it1e_end) {
188                                it1_index = it1.index1 ();
189                                it1e_index = it1e.index1 ();
190                            } else
191                                break;
192                        } else if (compare < 0) {
193                            increment (it1, it1_end, - compare);
194                            if (it1 != it1_end)
195                                it1_index = it1.index1 ();
196                            else
197                                break;
198                        } else if (compare > 0) {
199                            if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
200                                if (*it1e != value_type/*zero*/())
201                                    index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
202                            ++ it1e;
203                            if (it1e != it1e_end)
204                                it1e_index = it1e.index1 ();
205                            else
206                                break;
207                        }
208                    }
209                }
210                while (it1e != it1e_end) {
211                    if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
212                        if (*it1e != value_type/*zero*/())
213                            index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
214                    ++ it1e;
215                }
216                ++ it2, ++ it2e;
217            } else if (compare < 0) {
218                increment (it2, it2_end, - compare);
219            } else if (compare > 0) {
220#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
221                typename E::const_iterator1 it1e (it2e.begin ());
222                typename E::const_iterator1 it1e_end (it2e.end ());
223#else
224                typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
225                typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
226#endif
227                while (it1e != it1e_end) {
228                    if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
229                        if (*it1e != value_type/*zero*/())
230                            index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
231                    ++ it1e;
232                }
233                ++ it2e;
234            }
235        }
236        while (it2e != it2e_end) {
237#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
238            typename E::const_iterator1 it1e (it2e.begin ());
239            typename E::const_iterator1 it1e_end (it2e.end ());
240#else
241            typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
242            typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
243#endif
244            while (it1e != it1e_end) {
245                if (conformant_restrict_type::other (it1e.index1 (), it1e.index2 ()))
246                    if (*it1e != value_type/*zero*/())
247                        index.push_back (std::pair<size_type, size_type> (it1e.index1 (), it1e.index2 ()));
248                ++ it1e;
249            }
250            ++ it2e;
251        }
252        // ISSUE proxies require insert_element
253        for (size_type k = 0; k < index.size (); ++ k)
254            m (index [k].first, index [k].second) = value_type/*zero*/();
255    }
256
257}//namespace detail
258
259
260    // Explicitly iterating row major
261    template<template <class T1, class T2> class F, class M, class T>
262    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
263    void iterating_matrix_assign_scalar (M &m, const T &t, row_major_tag) {
264        typedef F<typename M::iterator2::reference, T> functor_type;
265        typedef typename M::difference_type difference_type;
266        difference_type size1 (m.size1 ());
267        difference_type size2 (m.size2 ());
268        typename M::iterator1 it1 (m.begin1 ());
269        BOOST_UBLAS_CHECK (size2 == 0 || m.end1 () - it1 == size1, bad_size ());
270        while (-- size1 >= 0) {
271#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
272            typename M::iterator2 it2 (it1.begin ());
273#else
274            typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
275#endif
276            BOOST_UBLAS_CHECK (it1.end () - it2 == size2, bad_size ());
277            difference_type temp_size2 (size2);
278#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
279            while (-- temp_size2 >= 0)
280                functor_type::apply (*it2, t), ++ it2;
281#else
282            DD (temp_size2, 4, r, (functor_type::apply (*it2, t), ++ it2));
283#endif
284            ++ it1;
285        }
286    }
287    // Explicitly iterating column major
288    template<template <class T1, class T2> class F, class M, class T>
289    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
290    void iterating_matrix_assign_scalar (M &m, const T &t, column_major_tag) {
291        typedef F<typename M::iterator1::reference, T> functor_type;
292        typedef typename M::difference_type difference_type;
293        difference_type size2 (m.size2 ());
294        difference_type size1 (m.size1 ());
295        typename M::iterator2 it2 (m.begin2 ());
296        BOOST_UBLAS_CHECK (size1 == 0 || m.end2 () - it2 == size2, bad_size ());
297        while (-- size2 >= 0) {
298#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
299            typename M::iterator1 it1 (it2.begin ());
300#else
301            typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
302#endif
303            BOOST_UBLAS_CHECK (it2.end () - it1 == size1, bad_size ());
304            difference_type temp_size1 (size1);
305#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
306            while (-- temp_size1 >= 0)
307                functor_type::apply (*it1, t), ++ it1;
308#else
309            DD (temp_size1, 4, r, (functor_type::apply (*it1, t), ++ it1));
310#endif
311            ++ it2;
312        }
313    }
314    // Explicitly indexing row major
315    template<template <class T1, class T2> class F, class M, class T>
316    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
317    void indexing_matrix_assign_scalar (M &m, const T &t, row_major_tag) {
318        typedef F<typename M::reference, T> functor_type;
319        typedef typename M::size_type size_type;
320        size_type size1 (m.size1 ());
321        size_type size2 (m.size2 ());
322        for (size_type i = 0; i < size1; ++ i) {
323#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
324            for (size_type j = 0; j < size2; ++ j)
325                functor_type::apply (m (i, j), t);
326#else
327            size_type j (0);
328            DD (size2, 4, r, (functor_type::apply (m (i, j), t), ++ j));
329#endif
330        }
331    }
332    // Explicitly indexing column major
333    template<template <class T1, class T2> class F, class M, class T>
334    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
335    void indexing_matrix_assign_scalar (M &m, const T &t, column_major_tag) {
336        typedef F<typename M::reference, T> functor_type;
337        typedef typename M::size_type size_type;
338        size_type size2 (m.size2 ());
339        size_type size1 (m.size1 ());
340        for (size_type j = 0; j < size2; ++ j) {
341#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
342            for (size_type i = 0; i < size1; ++ i)
343                functor_type::apply (m (i, j), t);
344#else
345            size_type i (0);
346            DD (size1, 4, r, (functor_type::apply (m (i, j), t), ++ i));
347#endif
348        }
349    }
350
351    // Dense (proxy) case
352    template<template <class T1, class T2> class F, class M, class T, class C>
353    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
354    void matrix_assign_scalar (M &m, const T &t, dense_proxy_tag, C) {
355        typedef C orientation_category;
356#ifdef BOOST_UBLAS_USE_INDEXING
357        indexing_matrix_assign_scalar<F> (m, t, orientation_category ());
358#elif BOOST_UBLAS_USE_ITERATING
359        iterating_matrix_assign_scalar<F> (m, t, orientation_category ());
360#else
361        typedef typename M::size_type size_type;
362        size_type size1 (m.size1 ());
363        size_type size2 (m.size2 ());
364        if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
365            size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
366            iterating_matrix_assign_scalar<F> (m, t, orientation_category ());
367        else
368            indexing_matrix_assign_scalar<F> (m, t, orientation_category ());
369#endif
370    }
371    // Packed (proxy) row major case
372    template<template <class T1, class T2> class F, class M, class T>
373    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
374    void matrix_assign_scalar (M &m, const T &t, packed_proxy_tag, row_major_tag) {
375        typedef F<typename M::iterator2::reference, T> functor_type;
376        typedef typename M::difference_type difference_type;
377        typename M::iterator1 it1 (m.begin1 ());
378        difference_type size1 (m.end1 () - it1);
379        while (-- size1 >= 0) {
380#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
381            typename M::iterator2 it2 (it1.begin ());
382            difference_type size2 (it1.end () - it2);
383#else
384            typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
385            difference_type size2 (end (it1, iterator1_tag ()) - it2);
386#endif
387            while (-- size2 >= 0)
388                functor_type::apply (*it2, t), ++ it2;
389            ++ it1;
390        }
391    }
392    // Packed (proxy) column major case
393    template<template <class T1, class T2> class F, class M, class T>
394    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
395    void matrix_assign_scalar (M &m, const T &t, packed_proxy_tag, column_major_tag) {
396        typedef F<typename M::iterator1::reference, T> functor_type;
397        typedef typename M::difference_type difference_type;
398        typename M::iterator2 it2 (m.begin2 ());
399        difference_type size2 (m.end2 () - it2);
400        while (-- size2 >= 0) {
401#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
402            typename M::iterator1 it1 (it2.begin ());
403            difference_type size1 (it2.end () - it1);
404#else
405            typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
406            difference_type size1 (end (it2, iterator2_tag ()) - it1);
407#endif
408            while (-- size1 >= 0)
409                functor_type::apply (*it1, t), ++ it1;
410            ++ it2;
411        }
412    }
413    // Sparse (proxy) row major case
414    template<template <class T1, class T2> class F, class M, class T>
415    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
416    void matrix_assign_scalar (M &m, const T &t, sparse_proxy_tag, row_major_tag) {
417        typedef F<typename M::iterator2::reference, T> functor_type;
418        typename M::iterator1 it1 (m.begin1 ());
419        typename M::iterator1 it1_end (m.end1 ());
420        while (it1 != it1_end) {
421#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
422            typename M::iterator2 it2 (it1.begin ());
423            typename M::iterator2 it2_end (it1.end ());
424#else
425            typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
426            typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
427#endif
428            while (it2 != it2_end)
429                functor_type::apply (*it2, t), ++ it2;
430            ++ it1;
431        }
432    }
433    // Sparse (proxy) column major case
434    template<template <class T1, class T2> class F, class M, class T>
435    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
436    void matrix_assign_scalar (M &m, const T &t, sparse_proxy_tag, column_major_tag) {
437        typedef F<typename M::iterator1::reference, T> functor_type;
438        typename M::iterator2 it2 (m.begin2 ());
439        typename M::iterator2 it2_end (m.end2 ());
440        while (it2 != it2_end) {
441#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
442            typename M::iterator1 it1 (it2.begin ());
443            typename M::iterator1 it1_end (it2.end ());
444#else
445            typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
446            typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
447#endif
448            while (it1 != it1_end)
449                functor_type::apply (*it1, t), ++ it1;
450            ++ it2;
451        }
452    }
453
454    // Dispatcher
455    template<template <class T1, class T2> class F, class M, class T>
456    BOOST_UBLAS_INLINE
457    void matrix_assign_scalar (M &m, const T &t) {
458        typedef typename M::storage_category storage_category;
459        typedef typename M::orientation_category orientation_category;
460        matrix_assign_scalar<F> (m, t, storage_category (), orientation_category ());
461    }
462
463    template<class SC, bool COMPUTED, class RI1, class RI2>
464    struct matrix_assign_traits {
465        typedef SC storage_category;
466    };
467
468    template<bool COMPUTED>
469    struct matrix_assign_traits<dense_tag, COMPUTED, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
470        typedef packed_tag storage_category;
471    };
472    template<>
473    struct matrix_assign_traits<dense_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
474        typedef sparse_tag storage_category;
475    };
476    template<>
477    struct matrix_assign_traits<dense_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
478        typedef sparse_proxy_tag storage_category;
479    };
480
481    template<bool COMPUTED>
482    struct matrix_assign_traits<dense_proxy_tag, COMPUTED, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
483        typedef packed_proxy_tag storage_category;
484    };
485    template<bool COMPUTED>
486    struct matrix_assign_traits<dense_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
487        typedef sparse_proxy_tag storage_category;
488    };
489
490    template<>
491    struct matrix_assign_traits<packed_tag, false, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
492        typedef sparse_tag storage_category;
493    };
494    template<>
495    struct matrix_assign_traits<packed_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
496        typedef sparse_proxy_tag storage_category;
497    };
498
499    template<bool COMPUTED>
500    struct matrix_assign_traits<packed_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
501        typedef sparse_proxy_tag storage_category;
502    };
503
504    template<>
505    struct matrix_assign_traits<sparse_tag, true, dense_random_access_iterator_tag, dense_random_access_iterator_tag> {
506        typedef sparse_proxy_tag storage_category;
507    };
508    template<>
509    struct matrix_assign_traits<sparse_tag, true, packed_random_access_iterator_tag, packed_random_access_iterator_tag> {
510        typedef sparse_proxy_tag storage_category;
511    };
512    template<>
513    struct matrix_assign_traits<sparse_tag, true, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
514        typedef sparse_proxy_tag storage_category;
515    };
516
517    // Explicitly iterating row major
518    template<template <class T1, class T2> class F, class M, class E>
519    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
520    void iterating_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
521        typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
522        typedef typename M::difference_type difference_type;
523        difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
524        difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
525        typename M::iterator1 it1 (m.begin1 ());
526        BOOST_UBLAS_CHECK (size2 == 0 || m.end1 () - it1 == size1, bad_size ());
527        typename E::const_iterator1 it1e (e ().begin1 ());
528        BOOST_UBLAS_CHECK (size2 == 0 || e ().end1 () - it1e == size1, bad_size ());
529        while (-- size1 >= 0) {
530#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
531            typename M::iterator2 it2 (it1.begin ());
532            typename E::const_iterator2 it2e (it1e.begin ());
533#else
534            typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
535            typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
536#endif
537            BOOST_UBLAS_CHECK (it1.end () - it2 == size2, bad_size ());
538            BOOST_UBLAS_CHECK (it1e.end () - it2e == size2, bad_size ());
539            difference_type temp_size2 (size2);
540#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
541            while (-- temp_size2 >= 0)
542                functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
543#else
544            DD (temp_size2, 2, r, (functor_type::apply (*it2, *it2e), ++ it2, ++ it2e));
545#endif
546            ++ it1, ++ it1e;
547        }
548    }
549    // Explicitly iterating column major
550    template<template <class T1, class T2> class F, class M, class E>
551    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
552    void iterating_matrix_assign (M &m, const matrix_expression<E> &e, column_major_tag) {
553        typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
554        typedef typename M::difference_type difference_type;
555        difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
556        difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
557        typename M::iterator2 it2 (m.begin2 ());
558        BOOST_UBLAS_CHECK (size1 == 0 || m.end2 () - it2 == size2, bad_size ());
559        typename E::const_iterator2 it2e (e ().begin2 ());
560        BOOST_UBLAS_CHECK (size1 == 0 || e ().end2 () - it2e == size2, bad_size ());
561        while (-- size2 >= 0) {
562#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
563            typename M::iterator1 it1 (it2.begin ());
564            typename E::const_iterator1 it1e (it2e.begin ());
565#else
566            typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
567            typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
568#endif
569            BOOST_UBLAS_CHECK (it2.end () - it1 == size1, bad_size ());
570            BOOST_UBLAS_CHECK (it2e.end () - it1e == size1, bad_size ());
571            difference_type temp_size1 (size1);
572#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
573            while (-- temp_size1 >= 0)
574                functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
575#else
576            DD (temp_size1, 2, r, (functor_type::apply (*it1, *it1e), ++ it1, ++ it1e));
577#endif
578            ++ it2, ++ it2e;
579        }
580    }
581    // Explicitly indexing row major
582    template<template <class T1, class T2> class F, class M, class E>
583    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
584    void indexing_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
585        typedef F<typename M::reference, typename E::value_type> functor_type;
586        typedef typename M::size_type size_type;
587        size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
588        size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
589        for (size_type i = 0; i < size1; ++ i) {
590#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
591            for (size_type j = 0; j < size2; ++ j)
592                functor_type::apply (m (i, j), e () (i, j));
593#else
594            size_type j (0);
595            DD (size2, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ j));
596#endif
597        }
598    }
599    // Explicitly indexing column major
600    template<template <class T1, class T2> class F, class M, class E>
601    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
602    void indexing_matrix_assign (M &m, const matrix_expression<E> &e, column_major_tag) {
603        typedef F<typename M::reference, typename E::value_type> functor_type;
604        typedef typename M::size_type size_type;
605        size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
606        size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
607        for (size_type j = 0; j < size2; ++ j) {
608#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
609            for (size_type i = 0; i < size1; ++ i)
610                functor_type::apply (m (i, j), e () (i, j));
611#else
612            size_type i (0);
613            DD (size1, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ i));
614#endif
615        }
616    }
617
618    // Dense (proxy) case
619    template<template <class T1, class T2> class F, class R, class M, class E, class C>
620    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
621    void matrix_assign (M &m, const matrix_expression<E> &e, dense_proxy_tag, C) {
622        // R unnecessary, make_conformant not required
623        typedef C orientation_category;
624#ifdef BOOST_UBLAS_USE_INDEXING
625        indexing_matrix_assign<F> (m, e, orientation_category ());
626#elif BOOST_UBLAS_USE_ITERATING
627        iterating_matrix_assign<F> (m, e, orientation_category ());
628#else
629        typedef typename M::difference_type difference_type;
630        size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
631        size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
632        if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
633            size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
634            iterating_matrix_assign<F> (m, e, orientation_category ());
635        else
636            indexing_matrix_assign<F> (m, e, orientation_category ());
637#endif
638    }
639    // Packed (proxy) row major case
640    template<template <class T1, class T2> class F, class R, class M, class E>
641    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
642    void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
643        typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
644        // R unnecessary, make_conformant not required
645        typedef typename M::difference_type difference_type;
646        typedef typename M::value_type value_type;
647        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
648        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
649#if BOOST_UBLAS_TYPE_CHECK
650        matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
651        indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ());
652        indexing_matrix_assign<F> (cm, e, row_major_tag ());
653#endif
654        typename M::iterator1 it1 (m.begin1 ());
655        typename M::iterator1 it1_end (m.end1 ());
656        typename E::const_iterator1 it1e (e ().begin1 ());
657        typename E::const_iterator1 it1e_end (e ().end1 ());
658        difference_type it1_size (it1_end - it1);
659        difference_type it1e_size (it1e_end - it1e);
660        difference_type diff1 (0);
661        if (it1_size > 0 && it1e_size > 0)
662            diff1 = it1.index1 () - it1e.index1 ();
663        if (diff1 != 0) {
664            difference_type size1 = (std::min) (diff1, it1e_size);
665            if (size1 > 0) {
666                it1e += size1;
667                it1e_size -= size1;
668                diff1 -= size1;
669            }
670            size1 = (std::min) (- diff1, it1_size);
671            if (size1 > 0) {
672                it1_size -= size1;
673                if (!functor_type::computed) {
674                    while (-- size1 >= 0) { // zeroing
675#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
676                        typename M::iterator2 it2 (it1.begin ());
677                        typename M::iterator2 it2_end (it1.end ());
678#else
679                        typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
680                        typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
681#endif
682                        difference_type size2 (it2_end - it2);
683                        while (-- size2 >= 0)
684                            functor_type::apply (*it2, value_type/*zero*/()), ++ it2;
685                        ++ it1;
686                    }
687                } else {
688                    it1 += size1;
689                }
690                diff1 += size1;
691            }
692        }
693        difference_type size1 ((std::min) (it1_size, it1e_size));
694        it1_size -= size1;
695        it1e_size -= size1;
696        while (-- size1 >= 0) {
697#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
698            typename M::iterator2 it2 (it1.begin ());
699            typename M::iterator2 it2_end (it1.end ());
700            typename E::const_iterator2 it2e (it1e.begin ());
701            typename E::const_iterator2 it2e_end (it1e.end ());
702#else
703            typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
704            typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
705            typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
706            typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
707#endif
708            difference_type it2_size (it2_end - it2);
709            difference_type it2e_size (it2e_end - it2e);
710            difference_type diff2 (0);
711            if (it2_size > 0 && it2e_size > 0) {
712                diff2 = it2.index2 () - it2e.index2 ();
713                difference_type size2 = (std::min) (diff2, it2e_size);
714                if (size2 > 0) {
715                    it2e += size2;
716                    it2e_size -= size2;
717                    diff2 -= size2;
718                }
719                size2 = (std::min) (- diff2, it2_size);
720                if (size2 > 0) {
721                    it2_size -= size2;
722                    if (!functor_type::computed) {
723                        while (-- size2 >= 0)   // zeroing
724                            functor_type::apply (*it2, value_type/*zero*/()), ++ it2;
725                    } else {
726                        it2 += size2;
727                    }
728                    diff2 += size2;
729                }
730            }
731            difference_type size2 ((std::min) (it2_size, it2e_size));
732            it2_size -= size2;
733            it2e_size -= size2;
734            while (-- size2 >= 0)
735                functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
736            size2 = it2_size;
737            if (!functor_type::computed) {
738                while (-- size2 >= 0)   // zeroing
739                    functor_type::apply (*it2, value_type/*zero*/()), ++ it2;
740            } else {
741                it2 += size2;
742            }
743            ++ it1, ++ it1e;
744        }
745        size1 = it1_size;
746        if (!functor_type::computed) {
747            while (-- size1 >= 0) { // zeroing
748#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
749                typename M::iterator2 it2 (it1.begin ());
750                typename M::iterator2 it2_end (it1.end ());
751#else
752                typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
753                typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
754#endif
755                difference_type size2 (it2_end - it2);
756                while (-- size2 >= 0)
757                    functor_type::apply (*it2, value_type/*zero*/()), ++ it2;
758                ++ it1;
759            }
760        } else {
761            it1 += size1;
762        }
763#if BOOST_UBLAS_TYPE_CHECK
764        if (! disable_type_check<bool>::value)
765            BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
766#endif
767    }
768    // Packed (proxy) column major case
769    template<template <class T1, class T2> class F, class R, class M, class E>
770    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
771    void matrix_assign (M &m, const matrix_expression<E> &e, packed_proxy_tag, column_major_tag) {
772        typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
773        // R unnecessary, make_conformant not required
774        typedef typename M::difference_type difference_type;
775        typedef typename M::value_type value_type;
776        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
777        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
778#if BOOST_UBLAS_TYPE_CHECK
779        matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
780        indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ());
781        indexing_matrix_assign<F> (cm, e, column_major_tag ());
782#endif
783        typename M::iterator2 it2 (m.begin2 ());
784        typename M::iterator2 it2_end (m.end2 ());
785        typename E::const_iterator2 it2e (e ().begin2 ());
786        typename E::const_iterator2 it2e_end (e ().end2 ());
787        difference_type it2_size (it2_end - it2);
788        difference_type it2e_size (it2e_end - it2e);
789        difference_type diff2 (0);
790        if (it2_size > 0 && it2e_size > 0)
791            diff2 = it2.index2 () - it2e.index2 ();
792        if (diff2 != 0) {
793            difference_type size2 = (std::min) (diff2, it2e_size);
794            if (size2 > 0) {
795                it2e += size2;
796                it2e_size -= size2;
797                diff2 -= size2;
798            }
799            size2 = (std::min) (- diff2, it2_size);
800            if (size2 > 0) {
801                it2_size -= size2;
802                if (!functor_type::computed) {
803                    while (-- size2 >= 0) { // zeroing
804#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
805                        typename M::iterator1 it1 (it2.begin ());
806                        typename M::iterator1 it1_end (it2.end ());
807#else
808                        typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
809                        typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
810#endif
811                        difference_type size1 (it1_end - it1);
812                        while (-- size1 >= 0)
813                            functor_type::apply (*it1, value_type/*zero*/()), ++ it1;
814                        ++ it2;
815                    }
816                } else {
817                    it2 += size2;
818                }
819                diff2 += size2;
820            }
821        }
822        difference_type size2 ((std::min) (it2_size, it2e_size));
823        it2_size -= size2;
824        it2e_size -= size2;
825        while (-- size2 >= 0) {
826#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
827            typename M::iterator1 it1 (it2.begin ());
828            typename M::iterator1 it1_end (it2.end ());
829            typename E::const_iterator1 it1e (it2e.begin ());
830            typename E::const_iterator1 it1e_end (it2e.end ());
831#else
832            typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
833            typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
834            typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
835            typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
836#endif
837            difference_type it1_size (it1_end - it1);
838            difference_type it1e_size (it1e_end - it1e);
839            difference_type diff1 (0);
840            if (it1_size > 0 && it1e_size > 0) {
841                diff1 = it1.index1 () - it1e.index1 ();
842                difference_type size1 = (std::min) (diff1, it1e_size);
843                if (size1 > 0) {
844                    it1e += size1;
845                    it1e_size -= size1;
846                    diff1 -= size1;
847                }
848                size1 = (std::min) (- diff1, it1_size);
849                if (size1 > 0) {
850                    it1_size -= size1;
851                    if (!functor_type::computed) {
852                        while (-- size1 >= 0)   // zeroing
853                            functor_type::apply (*it1, value_type/*zero*/()), ++ it1;
854                    } else {
855                        it1 += size1;
856                    }
857                    diff1 += size1;
858                }
859            }
860            difference_type size1 ((std::min) (it1_size, it1e_size));
861            it1_size -= size1;
862            it1e_size -= size1;
863            while (-- size1 >= 0)
864                functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
865            size1 = it1_size;
866            if (!functor_type::computed) {
867                while (-- size1 >= 0)   // zeroing
868                    functor_type::apply (*it1, value_type/*zero*/()), ++ it1;
869            } else {
870                it1 += size1;
871            }
872            ++ it2, ++ it2e;
873        }
874        size2 = it2_size;
875        if (!functor_type::computed) {
876            while (-- size2 >= 0) { // zeroing
877#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
878                typename M::iterator1 it1 (it2.begin ());
879                typename M::iterator1 it1_end (it2.end ());
880#else
881                typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
882                typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
883#endif
884                difference_type size1 (it1_end - it1);
885                while (-- size1 >= 0)
886                    functor_type::apply (*it1, value_type/*zero*/()), ++ it1;
887                ++ it2;
888            }
889        } else {
890            it2 += size2;
891        }
892#if BOOST_UBLAS_TYPE_CHECK
893        if (! disable_type_check<bool>::value)
894            BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
895#endif
896    }
897    // Sparse row major case
898    template<template <class T1, class T2> class F, class R, class M, class E>
899    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
900    void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, row_major_tag) {
901        typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
902        // R unnecessary, make_conformant not required
903        BOOST_STATIC_ASSERT ((!functor_type::computed));
904        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
905        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
906        typedef typename M::value_type value_type;
907        // Sparse type has no numeric constraints to check
908
909        m.clear ();
910        typename E::const_iterator1 it1e (e ().begin1 ());
911        typename E::const_iterator1 it1e_end (e ().end1 ());
912        while (it1e != it1e_end) {
913#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
914            typename E::const_iterator2 it2e (it1e.begin ());
915            typename E::const_iterator2 it2e_end (it1e.end ());
916#else
917            typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
918            typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
919#endif
920            while (it2e != it2e_end) {
921                value_type t (*it2e);
922                if (t != value_type/*zero*/())
923                    m.insert_element (it2e.index1 (), it2e.index2 (), t);
924                ++ it2e;
925            }
926            ++ it1e;
927        }
928    }
929    // Sparse column major case
930    template<template <class T1, class T2> class F, class R, class M, class E>
931    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
932    void matrix_assign (M &m, const matrix_expression<E> &e, sparse_tag, column_major_tag) {
933        typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
934        // R unnecessary, make_conformant not required
935        BOOST_STATIC_ASSERT ((!functor_type::computed));
936        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
937        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
938        typedef typename M::value_type value_type;
939        // Sparse type has no numeric constraints to check
940
941        m.clear ();
942        typename E::const_iterator2 it2e (e ().begin2 ());
943        typename E::const_iterator2 it2e_end (e ().end2 ());
944        while (it2e != it2e_end) {
945#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
946            typename E::const_iterator1 it1e (it2e.begin ());
947            typename E::const_iterator1 it1e_end (it2e.end ());
948#else
949            typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
950            typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
951#endif
952            while (it1e != it1e_end) {
953                value_type t (*it1e);
954                if (t != value_type/*zero*/())
955                    m.insert_element (it1e.index1 (), it1e.index2 (), t);
956                ++ it1e;
957            }
958            ++ it2e;
959        }
960    }
961    // Sparse proxy or functional row major case
962    template<template <class T1, class T2> class F, class R, class M, class E>
963    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
964    void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) {
965        typedef F<typename M::iterator2::reference, typename E::value_type> functor_type;
966        typedef R conformant_restrict_type;
967        typedef typename M::size_type size_type;
968        typedef typename M::difference_type difference_type;
969        typedef typename M::value_type value_type;
970        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
971        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
972#if BOOST_UBLAS_TYPE_CHECK
973        matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
974        indexing_matrix_assign<scalar_assign> (cm, m, row_major_tag ());
975        indexing_matrix_assign<F> (cm, e, row_major_tag ());
976#endif
977        detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
978
979        typename M::iterator1 it1 (m.begin1 ());
980        typename M::iterator1 it1_end (m.end1 ());
981        typename E::const_iterator1 it1e (e ().begin1 ());
982        typename E::const_iterator1 it1e_end (e ().end1 ());
983        while (it1 != it1_end && it1e != it1e_end) {
984            difference_type compare = it1.index1 () - it1e.index1 ();
985            if (compare == 0) {
986#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
987                typename M::iterator2 it2 (it1.begin ());
988                typename M::iterator2 it2_end (it1.end ());
989                typename E::const_iterator2 it2e (it1e.begin ());
990                typename E::const_iterator2 it2e_end (it1e.end ());
991#else
992                typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
993                typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
994                typename E::const_iterator2 it2e (begin (it1e, iterator1_tag ()));
995                typename E::const_iterator2 it2e_end (end (it1e, iterator1_tag ()));
996#endif
997                if (it2 != it2_end && it2e != it2e_end) {
998                    size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
999                    while (true) {
1000                        difference_type compare = it2_index - it2e_index;
1001                        if (compare == 0) {
1002                            functor_type::apply (*it2, *it2e);
1003                            ++ it2, ++ it2e;
1004                            if (it2 != it2_end && it2e != it2e_end) {
1005                                it2_index = it2.index2 ();
1006                                it2e_index = it2e.index2 ();
1007                            } else
1008                                break;
1009                        } else if (compare < 0) {
1010                            if (!functor_type::computed) {
1011                                functor_type::apply (*it2, value_type/*zero*/());
1012                                ++ it2;
1013                            } else
1014                                increment (it2, it2_end, - compare);
1015                            if (it2 != it2_end)
1016                                it2_index = it2.index2 ();
1017                            else
1018                                break;
1019                        } else if (compare > 0) {
1020                            increment (it2e, it2e_end, compare);
1021                            if (it2e != it2e_end)
1022                                it2e_index = it2e.index2 ();
1023                            else
1024                                break;
1025                        }
1026                    }
1027                }
1028                if (!functor_type::computed) {
1029                    while (it2 != it2_end) {    // zeroing
1030                        functor_type::apply (*it2, value_type/*zero*/());
1031                        ++ it2;
1032                    }
1033                } else {
1034                    it2 = it2_end;
1035                }
1036                ++ it1, ++ it1e;
1037            } else if (compare < 0) {
1038                if (!functor_type::computed) {
1039#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1040                    typename M::iterator2 it2 (it1.begin ());
1041                    typename M::iterator2 it2_end (it1.end ());
1042#else
1043                    typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1044                    typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1045#endif
1046                    while (it2 != it2_end) {    // zeroing
1047                        functor_type::apply (*it2, value_type/*zero*/());
1048                        ++ it2;
1049                    }
1050                    ++ it1;
1051                } else {
1052                    increment (it1, it1_end, - compare);
1053                }
1054            } else if (compare > 0) {
1055                increment (it1e, it1e_end, compare);
1056            }
1057        }
1058        if (!functor_type::computed) {
1059            while (it1 != it1_end) {
1060#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1061                typename M::iterator2 it2 (it1.begin ());
1062                typename M::iterator2 it2_end (it1.end ());
1063#else
1064                typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1065                typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1066#endif
1067                while (it2 != it2_end) {    // zeroing
1068                    functor_type::apply (*it2, value_type/*zero*/());
1069                    ++ it2;
1070                }
1071                ++ it1;
1072            }
1073        } else {
1074            it1 = it1_end;
1075        }
1076#if BOOST_UBLAS_TYPE_CHECK
1077        if (! disable_type_check<bool>::value)
1078            BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
1079#endif
1080    }
1081    // Sparse proxy or functional column major case
1082    template<template <class T1, class T2> class F, class R, class M, class E>
1083    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1084    void matrix_assign (M &m, const matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) {
1085        typedef F<typename M::iterator1::reference, typename E::value_type> functor_type;
1086        typedef R conformant_restrict_type;
1087        typedef typename M::size_type size_type;
1088        typedef typename M::difference_type difference_type;
1089        typedef typename M::value_type value_type;
1090        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
1091        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
1092#if BOOST_UBLAS_TYPE_CHECK
1093        matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
1094        indexing_matrix_assign<scalar_assign> (cm, m, column_major_tag ());
1095        indexing_matrix_assign<F> (cm, e, column_major_tag ());
1096#endif
1097        detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
1098
1099        typename M::iterator2 it2 (m.begin2 ());
1100        typename M::iterator2 it2_end (m.end2 ());
1101        typename E::const_iterator2 it2e (e ().begin2 ());
1102        typename E::const_iterator2 it2e_end (e ().end2 ());
1103        while (it2 != it2_end && it2e != it2e_end) {
1104            difference_type compare = it2.index2 () - it2e.index2 ();
1105            if (compare == 0) {
1106#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1107                typename M::iterator1 it1 (it2.begin ());
1108                typename M::iterator1 it1_end (it2.end ());
1109                typename E::const_iterator1 it1e (it2e.begin ());
1110                typename E::const_iterator1 it1e_end (it2e.end ());
1111#else
1112                typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1113                typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1114                typename E::const_iterator1 it1e (begin (it2e, iterator2_tag ()));
1115                typename E::const_iterator1 it1e_end (end (it2e, iterator2_tag ()));
1116#endif
1117                if (it1 != it1_end && it1e != it1e_end) {
1118                    size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
1119                    while (true) {
1120                        difference_type compare = it1_index - it1e_index;
1121                        if (compare == 0) {
1122                            functor_type::apply (*it1, *it1e);
1123                            ++ it1, ++ it1e;
1124                            if (it1 != it1_end && it1e != it1e_end) {
1125                                it1_index = it1.index1 ();
1126                                it1e_index = it1e.index1 ();
1127                            } else
1128                                break;
1129                        } else if (compare < 0) {
1130                            if (!functor_type::computed) {
1131                                functor_type::apply (*it1, value_type/*zero*/()); // zeroing
1132                                ++ it1;
1133                            } else
1134                                increment (it1, it1_end, - compare);
1135                            if (it1 != it1_end)
1136                                it1_index = it1.index1 ();
1137                            else
1138                                break;
1139                        } else if (compare > 0) {
1140                            increment (it1e, it1e_end, compare);
1141                            if (it1e != it1e_end)
1142                                it1e_index = it1e.index1 ();
1143                            else
1144                                break;
1145                        }
1146                    }
1147                }
1148                if (!functor_type::computed) {
1149                    while (it1 != it1_end) {    // zeroing
1150                        functor_type::apply (*it1, value_type/*zero*/());
1151                        ++ it1;
1152                    }
1153                } else {
1154                    it1 = it1_end;
1155                }
1156                ++ it2, ++ it2e;
1157            } else if (compare < 0) {
1158                if (!functor_type::computed) {
1159#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1160                    typename M::iterator1 it1 (it2.begin ());
1161                    typename M::iterator1 it1_end (it2.end ());
1162#else
1163                    typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1164                    typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1165#endif
1166                    while (it1 != it1_end) {    // zeroing
1167                        functor_type::apply (*it1, value_type/*zero*/());
1168                        ++ it1;
1169                    }
1170                    ++ it2;
1171                } else {
1172                    increment (it2, it2_end, - compare);
1173                }
1174            } else if (compare > 0) {
1175                increment (it2e, it2e_end, compare);
1176            }
1177        }
1178        if (!functor_type::computed) {
1179            while (it2 != it2_end) {
1180#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1181                typename M::iterator1 it1 (it2.begin ());
1182                typename M::iterator1 it1_end (it2.end ());
1183#else
1184                typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1185                typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1186#endif
1187                while (it1 != it1_end) {    // zeroing
1188                    functor_type::apply (*it1, value_type/*zero*/());
1189                    ++ it1;
1190                }
1191                ++ it2;
1192            }
1193        } else {
1194            it2 = it2_end;
1195        }
1196#if BOOST_UBLAS_TYPE_CHECK
1197        if (! disable_type_check<bool>::value)
1198            BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm), external_logic ());
1199#endif
1200    }
1201
1202    // Dispatcher
1203    template<template <class T1, class T2> class F, class M, class E>
1204    BOOST_UBLAS_INLINE
1205    void matrix_assign (M &m, const matrix_expression<E> &e) {
1206        typedef typename matrix_assign_traits<typename M::storage_category,
1207                                              F<typename M::reference, typename E::value_type>::computed,
1208                                              typename E::const_iterator1::iterator_category,
1209                                              typename E::const_iterator2::iterator_category>::storage_category storage_category;
1210        // give preference to matrix M's orientation if known
1211        typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
1212                                          typename E::orientation_category ,
1213                                          typename M::orientation_category >::type orientation_category;
1214        typedef basic_full<typename M::size_type> unrestricted;
1215        matrix_assign<F, unrestricted> (m, e, storage_category (), orientation_category ());
1216    }
1217    template<template <class T1, class T2> class F, class R, class M, class E>
1218    BOOST_UBLAS_INLINE
1219    void matrix_assign (M &m, const matrix_expression<E> &e) {
1220        typedef R conformant_restrict_type;
1221        typedef typename matrix_assign_traits<typename M::storage_category,
1222                                              F<typename M::reference, typename E::value_type>::computed,
1223                                              typename E::const_iterator1::iterator_category,
1224                                              typename E::const_iterator2::iterator_category>::storage_category storage_category;
1225        // give preference to matrix M's orientation if known
1226        typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
1227                                          typename E::orientation_category ,
1228                                          typename M::orientation_category >::type orientation_category;
1229        matrix_assign<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ());
1230    }
1231
1232    template<class SC, class RI1, class RI2>
1233    struct matrix_swap_traits {
1234        typedef SC storage_category;
1235    };
1236
1237    template<>
1238    struct matrix_swap_traits<dense_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
1239        typedef sparse_proxy_tag storage_category;
1240    };
1241
1242    template<>
1243    struct matrix_swap_traits<packed_proxy_tag, sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag> {
1244        typedef sparse_proxy_tag storage_category;
1245    };
1246
1247    // Dense (proxy) row major case
1248    template<template <class T1, class T2> class F, class R, class M, class E>
1249    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1250    void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, row_major_tag) {
1251        typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
1252        // R unnecessary, make_conformant not required
1253        typedef typename M::size_type size_type;
1254        typedef typename M::difference_type difference_type;
1255        typename M::iterator1 it1 (m.begin1 ());
1256        typename E::iterator1 it1e (e ().begin1 ());
1257        difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (e ().end1 () - it1e)));
1258        while (-- size1 >= 0) {
1259#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1260            typename M::iterator2 it2 (it1.begin ());
1261            typename E::iterator2 it2e (it1e.begin ());
1262            difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (it1e.end () - it2e)));
1263#else
1264            typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1265            typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1266            difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (end (it1e, iterator1_tag ()) - it2e)));
1267#endif
1268            while (-- size2 >= 0)
1269                functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
1270            ++ it1, ++ it1e;
1271        }
1272    }
1273    // Dense (proxy) column major case
1274    template<template <class T1, class T2> class F, class R, class M, class E>
1275    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1276    void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, column_major_tag) {
1277        typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
1278        // R unnecessary, make_conformant not required
1279        typedef typename M::size_type size_type;
1280        typedef typename M::difference_type difference_type;
1281        typename M::iterator2 it2 (m.begin2 ());
1282        typename E::iterator2 it2e (e ().begin2 ());
1283        difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (e ().end2 () - it2e)));
1284        while (-- size2 >= 0) {
1285#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1286            typename M::iterator1 it1 (it2.begin ());
1287            typename E::iterator1 it1e (it2e.begin ());
1288            difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (it2e.end () - it1e)));
1289#else
1290            typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1291            typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1292            difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (end (it2e, iterator2_tag ()) - it1e)));
1293#endif
1294            while (-- size1 >= 0)
1295                functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
1296            ++ it2, ++ it2e;
1297        }
1298    }
1299    // Packed (proxy) row major case
1300    template<template <class T1, class T2> class F, class R, class M, class E>
1301    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1302    void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
1303        typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
1304        // R unnecessary, make_conformant not required
1305        typedef typename M::size_type size_type;
1306        typedef typename M::difference_type difference_type;
1307        typename M::iterator1 it1 (m.begin1 ());
1308        typename E::iterator1 it1e (e ().begin1 ());
1309        difference_type size1 (BOOST_UBLAS_SAME (m.end1 () - it1, e ().end1 () - it1e));
1310        while (-- size1 >= 0) {
1311#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1312            typename M::iterator2 it2 (it1.begin ());
1313            typename E::iterator2 it2e (it1e.begin ());
1314            difference_type size2 (BOOST_UBLAS_SAME (it1.end () - it2, it1e.end () - it2e));
1315#else
1316            typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1317            typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1318            difference_type size2 (BOOST_UBLAS_SAME (end (it1, iterator1_tag ()) - it2, end (it1e, iterator1_tag ()) - it2e));
1319#endif
1320            while (-- size2 >= 0)
1321                functor_type::apply (*it2, *it2e), ++ it2, ++ it2e;
1322            ++ it1, ++ it1e;
1323        }
1324    }
1325    // Packed (proxy) column major case
1326    template<template <class T1, class T2> class F, class R, class M, class E>
1327    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1328    void matrix_swap (M &m, matrix_expression<E> &e, packed_proxy_tag, column_major_tag) {
1329        typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
1330        // R unnecessary, make_conformant not required
1331        typedef typename M::size_type size_type;
1332        typedef typename M::difference_type difference_type;
1333        typename M::iterator2 it2 (m.begin2 ());
1334        typename E::iterator2 it2e (e ().begin2 ());
1335        difference_type size2 (BOOST_UBLAS_SAME (m.end2 () - it2, e ().end2 () - it2e));
1336        while (-- size2 >= 0) {
1337#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1338            typename M::iterator1 it1 (it2.begin ());
1339            typename E::iterator1 it1e (it2e.begin ());
1340            difference_type size1 (BOOST_UBLAS_SAME (it2.end () - it1, it2e.end () - it1e));
1341#else
1342            typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1343            typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1344            difference_type size1 (BOOST_UBLAS_SAME (end (it2, iterator2_tag ()) - it1, end (it2e, iterator2_tag ()) - it1e));
1345#endif
1346            while (-- size1 >= 0)
1347                functor_type::apply (*it1, *it1e), ++ it1, ++ it1e;
1348            ++ it2, ++ it2e;
1349        }
1350    }
1351    // Sparse (proxy) row major case
1352    template<template <class T1, class T2> class F, class R, class M, class E>
1353    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1354    void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, row_major_tag) {
1355        typedef F<typename M::iterator2::reference, typename E::reference> functor_type;
1356        typedef R conformant_restrict_type;
1357        typedef typename M::size_type size_type;
1358        typedef typename M::difference_type difference_type;
1359        typedef typename M::value_type value_type;
1360        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
1361        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
1362
1363        detail::make_conformant (m, e, row_major_tag (), conformant_restrict_type ());
1364        // FIXME should be a seperate restriction for E
1365        detail::make_conformant (e (), m, row_major_tag (), conformant_restrict_type ());
1366
1367        typename M::iterator1 it1 (m.begin1 ());
1368        typename M::iterator1 it1_end (m.end1 ());
1369        typename E::iterator1 it1e (e ().begin1 ());
1370        typename E::iterator1 it1e_end (e ().end1 ());
1371        while (it1 != it1_end && it1e != it1e_end) {
1372            difference_type compare = it1.index1 () - it1e.index1 ();
1373            if (compare == 0) {
1374#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1375                typename M::iterator2 it2 (it1.begin ());
1376                typename M::iterator2 it2_end (it1.end ());
1377                typename E::iterator2 it2e (it1e.begin ());
1378                typename E::iterator2 it2e_end (it1e.end ());
1379#else
1380                typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1381                typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1382                typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1383                typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
1384#endif
1385                if (it2 != it2_end && it2e != it2e_end) {
1386                    size_type it2_index = it2.index2 (), it2e_index = it2e.index2 ();
1387                    while (true) {
1388                        difference_type compare = it2_index - it2e_index;
1389                        if (compare == 0) {
1390                            functor_type::apply (*it2, *it2e);
1391                            ++ it2, ++ it2e;
1392                            if (it2 != it2_end && it2e != it2e_end) {
1393                                it2_index = it2.index2 ();
1394                                it2e_index = it2e.index2 ();
1395                            } else
1396                                break;
1397                        } else if (compare < 0) {
1398                            increment (it2, it2_end, - compare);
1399                            if (it2 != it2_end)
1400                                it2_index = it2.index2 ();
1401                            else
1402                                break;
1403                        } else if (compare > 0) {
1404                            increment (it2e, it2e_end, compare);
1405                            if (it2e != it2e_end)
1406                                it2e_index = it2e.index2 ();
1407                            else
1408                                break;
1409                        }
1410                    }
1411                }
1412#if BOOST_UBLAS_TYPE_CHECK
1413                increment (it2e, it2e_end);
1414                increment (it2, it2_end);
1415#endif
1416                ++ it1, ++ it1e;
1417            } else if (compare < 0) {
1418#if BOOST_UBLAS_TYPE_CHECK
1419                while (it1.index1 () < it1e.index1 ()) {
1420#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1421                    typename M::iterator2 it2 (it1.begin ());
1422                    typename M::iterator2 it2_end (it1.end ());
1423#else
1424                    typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1425                    typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1426#endif
1427                    increment (it2, it2_end);
1428                    ++ it1;
1429                }
1430#else
1431                increment (it1, it1_end, - compare);
1432#endif
1433            } else if (compare > 0) {
1434#if BOOST_UBLAS_TYPE_CHECK
1435                while (it1e.index1 () < it1.index1 ()) {
1436#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1437                    typename E::iterator2 it2e (it1e.begin ());
1438                    typename E::iterator2 it2e_end (it1e.end ());
1439#else
1440                    typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1441                    typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
1442#endif
1443                    increment (it2e, it2e_end);
1444                    ++ it1e;
1445                }
1446#else
1447                increment (it1e, it1e_end, compare);
1448#endif
1449            }
1450        }
1451#if BOOST_UBLAS_TYPE_CHECK
1452        while (it1e != it1e_end) {
1453#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1454            typename E::iterator2 it2e (it1e.begin ());
1455            typename E::iterator2 it2e_end (it1e.end ());
1456#else
1457            typename E::iterator2 it2e (begin (it1e, iterator1_tag ()));
1458            typename E::iterator2 it2e_end (end (it1e, iterator1_tag ()));
1459#endif
1460            increment (it2e, it2e_end);
1461            ++ it1e;
1462        }
1463        while (it1 != it1_end) {
1464#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1465            typename M::iterator2 it2 (it1.begin ());
1466            typename M::iterator2 it2_end (it1.end ());
1467#else
1468            typename M::iterator2 it2 (begin (it1, iterator1_tag ()));
1469            typename M::iterator2 it2_end (end (it1, iterator1_tag ()));
1470#endif
1471            increment (it2, it2_end);
1472            ++ it1;
1473        }
1474#endif
1475    }
1476    // Sparse (proxy) column major case
1477    template<template <class T1, class T2> class F, class R, class M, class E>
1478    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1479    void matrix_swap (M &m, matrix_expression<E> &e, sparse_proxy_tag, column_major_tag) {
1480        typedef F<typename M::iterator1::reference, typename E::reference> functor_type;
1481        typedef R conformant_restrict_type;
1482        typedef typename M::size_type size_type;
1483        typedef typename M::difference_type difference_type;
1484        typedef typename M::value_type value_type;
1485        BOOST_UBLAS_CHECK (m.size1 () == e ().size1 (), bad_size ());
1486        BOOST_UBLAS_CHECK (m.size2 () == e ().size2 (), bad_size ());
1487
1488        detail::make_conformant (m, e, column_major_tag (), conformant_restrict_type ());
1489        // FIXME should be a seperate restriction for E
1490        detail::make_conformant (e (), m, column_major_tag (), conformant_restrict_type ());
1491
1492        typename M::iterator2 it2 (m.begin2 ());
1493        typename M::iterator2 it2_end (m.end2 ());
1494        typename E::iterator2 it2e (e ().begin2 ());
1495        typename E::iterator2 it2e_end (e ().end2 ());
1496        while (it2 != it2_end && it2e != it2e_end) {
1497            difference_type compare = it2.index2 () - it2e.index2 ();
1498            if (compare == 0) {
1499#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1500                typename M::iterator1 it1 (it2.begin ());
1501                typename M::iterator1 it1_end (it2.end ());
1502                typename E::iterator1 it1e (it2e.begin ());
1503                typename E::iterator1 it1e_end (it2e.end ());
1504#else
1505                typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1506                typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1507                typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1508                typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
1509#endif
1510                if (it1 != it1_end && it1e != it1e_end) {
1511                    size_type it1_index = it1.index1 (), it1e_index = it1e.index1 ();
1512                    while (true) {
1513                        difference_type compare = it1_index - it1e_index;
1514                        if (compare == 0) {
1515                            functor_type::apply (*it1, *it1e);
1516                            ++ it1, ++ it1e;
1517                            if (it1 != it1_end && it1e != it1e_end) {
1518                                it1_index = it1.index1 ();
1519                                it1e_index = it1e.index1 ();
1520                            } else
1521                                break;
1522                        }  else if (compare < 0) {
1523                            increment (it1, it1_end, - compare);
1524                            if (it1 != it1_end)
1525                                it1_index = it1.index1 ();
1526                            else
1527                                break;
1528                        } else if (compare > 0) {
1529                            increment (it1e, it1e_end, compare);
1530                            if (it1e != it1e_end)
1531                                it1e_index = it1e.index1 ();
1532                            else
1533                                break;
1534                        }
1535                    }
1536                }
1537#if BOOST_UBLAS_TYPE_CHECK
1538                increment (it1e, it1e_end);
1539                increment (it1, it1_end);
1540#endif
1541                ++ it2, ++ it2e;
1542            } else if (compare < 0) {
1543#if BOOST_UBLAS_TYPE_CHECK
1544                while (it2.index2 () < it2e.index2 ()) {
1545#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1546                    typename M::iterator1 it1 (it2.begin ());
1547                    typename M::iterator1 it1_end (it2.end ());
1548#else
1549                    typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1550                    typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1551#endif
1552                    increment (it1, it1_end);
1553                    ++ it2;
1554                }
1555#else
1556                increment (it2, it2_end, - compare);
1557#endif
1558            } else if (compare > 0) {
1559#if BOOST_UBLAS_TYPE_CHECK
1560                while (it2e.index2 () < it2.index2 ()) {
1561#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1562                    typename E::iterator1 it1e (it2e.begin ());
1563                    typename E::iterator1 it1e_end (it2e.end ());
1564#else
1565                    typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1566                    typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
1567#endif
1568                    increment (it1e, it1e_end);
1569                    ++ it2e;
1570                }
1571#else
1572                increment (it2e, it2e_end, compare);
1573#endif
1574            }
1575        }
1576#if BOOST_UBLAS_TYPE_CHECK
1577        while (it2e != it2e_end) {
1578#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1579            typename E::iterator1 it1e (it2e.begin ());
1580            typename E::iterator1 it1e_end (it2e.end ());
1581#else
1582            typename E::iterator1 it1e (begin (it2e, iterator2_tag ()));
1583            typename E::iterator1 it1e_end (end (it2e, iterator2_tag ()));
1584#endif
1585            increment (it1e, it1e_end);
1586            ++ it2e;
1587        }
1588        while (it2 != it2_end) {
1589#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1590            typename M::iterator1 it1 (it2.begin ());
1591            typename M::iterator1 it1_end (it2.end ());
1592#else
1593            typename M::iterator1 it1 (begin (it2, iterator2_tag ()));
1594            typename M::iterator1 it1_end (end (it2, iterator2_tag ()));
1595#endif
1596            increment (it1, it1_end);
1597            ++ it2;
1598        }
1599#endif
1600    }
1601
1602    // Dispatcher
1603    template<template <class T1, class T2> class F, class M, class E>
1604    BOOST_UBLAS_INLINE
1605    void matrix_swap (M &m, matrix_expression<E> &e) {
1606        typedef typename matrix_swap_traits<typename M::storage_category,
1607                                            typename E::const_iterator1::iterator_category,
1608                                            typename E::const_iterator2::iterator_category>::storage_category storage_category;
1609        // give preference to matrix M's orientation if known
1610        typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
1611                                          typename E::orientation_category ,
1612                                          typename M::orientation_category >::type orientation_category;
1613        typedef basic_full<typename M::size_type> unrestricted;
1614        matrix_swap<F, unrestricted> (m, e, storage_category (), orientation_category ());
1615    }
1616    template<template <class T1, class T2> class F, class R, class M, class E>
1617    BOOST_UBLAS_INLINE
1618    void matrix_swap (M &m, matrix_expression<E> &e) {
1619        typedef R conformant_restrict_type;
1620        typedef typename matrix_swap_traits<typename M::storage_category,
1621                                            typename E::const_iterator1::iterator_category,
1622                                            typename E::const_iterator2::iterator_category>::storage_category storage_category;
1623        // give preference to matrix M's orientation if known
1624        typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
1625                                          typename E::orientation_category ,
1626                                          typename M::orientation_category >::type orientation_category;
1627        matrix_swap<F, conformant_restrict_type> (m, e, storage_category (), orientation_category ());
1628    }
1629
1630}}}
1631
1632#endif
Note: See TracBrowser for help on using the repository browser.