← Back to branch summary

~twpol/dcplusplus/trunk

« back to all changes in this revision

Viewing changes to boost/boost/numeric/ublas/detail/matrix_assign.hpp

  • Committer: James Ross
  • Date: 2010-07-05 00:03:18 UTC
  • mfrom: (1524.1.650 dcplusplus)
  • Revision ID: silver@warwickcompsoc.co.uk-20100705000318-awwqm8ocpp5m47yz
Merged to trunk.

Show diffs side-by-side

added added

removed removed

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