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

Revision 857, 24.1 KB checked in by igarcia, 19 years ago (diff)
Line 
1//
2//  Copyright (c) 2000-2002
3//  Joerg Walter, Mathias Koch
4//
5//  Permission to use, copy, modify, distribute and sell this software
6//  and its documentation for any purpose is hereby granted without fee,
7//  provided that the above copyright notice appear in all copies and
8//  that both that copyright notice and this permission notice appear
9//  in supporting documentation.  The authors make no representations
10//  about the suitability of this software for any purpose.
11//  It is provided "as is" without express or implied warranty.
12//
13//  The authors gratefully acknowledge the support of
14//  GeNeSys mbH & Co. KG in producing this work.
15//
16
17#ifndef _BOOST_UBLAS_VECTOR_ASSIGN_
18#define _BOOST_UBLAS_VECTOR_ASSIGN_
19
20// Required for make_conformant storage
21#include <vector>
22
23// Iterators based on ideas of Jeremy Siek
24
25namespace boost { namespace numeric { namespace ublas {
26namespace detail {
27
28    // Weak equality check - useful to compare equality two arbitary vector expression results.
29    // Since the actual expressions are unknown, we check for and arbitary error bound
30    // on the relative error.
31    // For a linear expression the infinity norm makes sense as we do not know how the elements will be
32    // combined in the expression. False positive results are inevitable for arbirary expressions!
33    template<class E1, class E2, class S>
34    BOOST_UBLAS_INLINE
35    bool equals (const vector_expression<E1> &e1, const vector_expression<E2> &e2, S epsilon, S min_norm) {
36        return norm_inf (e1 - e2) < epsilon *
37               std::max<S> (std::max<S> (norm_inf (e1), norm_inf (e2)), min_norm);
38    }
39
40    template<class E1, class E2>
41    BOOST_UBLAS_INLINE
42    bool expression_type_check (const vector_expression<E1> &e1, const vector_expression<E2> &e2) {
43        typedef typename type_traits<typename promote_traits<typename E1::value_type,
44                                     typename E2::value_type>::promote_type>::real_type real_type;
45        return equals (e1, e2, BOOST_UBLAS_TYPE_CHECK_EPSILON, BOOST_UBLAS_TYPE_CHECK_MIN);
46    }
47
48
49    // Make sparse proxies conformant
50    template<class V, class E>
51    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
52    void make_conformant (V &v, const vector_expression<E> &e) {
53        BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
54        typedef typename V::size_type size_type;
55        typedef typename V::difference_type difference_type;
56        typedef typename V::value_type value_type;
57        // FIXME unbounded_array with push_back maybe better
58        std::vector<size_type> index;
59        typename V::iterator it (v.begin ());
60        typename V::iterator it_end (v.end ());
61        typename E::const_iterator ite (e ().begin ());
62        typename E::const_iterator ite_end (e ().end ());
63        if (it != it_end && ite != ite_end) {
64            size_type it_index = it.index (), ite_index = ite.index ();
65            while (true) {
66                difference_type compare = it_index - ite_index;
67                if (compare == 0) {
68                    ++ it, ++ ite;
69                    if (it != it_end && ite != ite_end) {
70                        it_index = it.index ();
71                        ite_index = ite.index ();
72                    } else
73                        break;
74                } else if (compare < 0) {
75                    increment (it, it_end, - compare);
76                    if (it != it_end)
77                        it_index = it.index ();
78                    else
79                        break;
80                } else if (compare > 0) {
81                    if (*ite != value_type/*zero*/())
82                        index.push_back (ite.index ());
83                    ++ ite;
84                    if (ite != ite_end)
85                        ite_index = ite.index ();
86                    else
87                        break;
88                }
89            }
90        }
91
92        while (ite != ite_end) {
93            if (*ite != value_type/*zero*/())
94                index.push_back (ite.index ());
95            ++ ite;
96        }
97        for (size_type k = 0; k < index.size (); ++ k)
98            v (index [k]) = value_type/*zero*/();
99    }
100
101}//namespace detail
102
103
104    // Explicitly iterating
105    template<template <class T1, class T2> class F, class V, class T>
106    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
107    void iterating_vector_assign_scalar (V &v, const T &t) {
108        typedef F<typename V::iterator::reference, T> functor_type;
109        typedef typename V::difference_type difference_type;
110        difference_type size (v.size ());
111        typename V::iterator it (v.begin ());
112        BOOST_UBLAS_CHECK (v.end () - it == size, bad_size ());
113#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
114        while (-- size >= 0)
115            functor_type::apply (*it, t), ++ it;
116#else
117        DD (size, 4, r, (functor_type::apply (*it, t), ++ it));
118#endif
119    }
120    // Explicitly case
121    template<template <class T1, class T2> class F, class V, class T>
122    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
123    void indexing_vector_assign_scalar (V &v, const T &t) {
124        typedef F<typename V::reference, T> functor_type;
125        typedef typename V::size_type size_type;
126        size_type size (v.size ());
127#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
128        for (size_type i = 0; i < size; ++ i)
129            functor_type::apply (v (i), t);
130#else
131        size_type i (0);
132        DD (size, 4, r, (functor_type::apply (v (i), t), ++ i));
133#endif
134    }
135
136    // Dense (proxy) case
137    template<template <class T1, class T2> class F, class V, class T>
138    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
139    void vector_assign_scalar (V &v, const T &t, dense_proxy_tag) {
140#ifdef BOOST_UBLAS_USE_INDEXING
141        indexing_vector_assign_scalar<F> (v, t);
142#elif BOOST_UBLAS_USE_ITERATING
143        iterating_vector_assign_scalar<F> (v, t);
144#else
145        typedef typename V::size_type size_type;
146        size_type size (v.size ());
147        if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
148            iterating_vector_assign_scalar<F> (v, t);
149        else
150            indexing_vector_assign_scalar<F> (v, t);
151#endif
152    }
153    // Packed (proxy) case
154    template<template <class T1, class T2> class F, class V, class T>
155    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
156    void vector_assign_scalar (V &v, const T &t, packed_proxy_tag) {
157        typedef F<typename V::iterator::reference, T> functor_type;
158        typedef typename V::difference_type difference_type;
159        typename V::iterator it (v.begin ());
160        difference_type size (v.end () - it);
161        while (-- size >= 0)
162            functor_type::apply (*it, t), ++ it;
163    }
164    // Sparse (proxy) case
165    template<template <class T1, class T2> class F, class V, class T>
166    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
167    void vector_assign_scalar (V &v, const T &t, sparse_proxy_tag) {
168        typedef F<typename V::iterator::reference, T> functor_type;
169        typename V::iterator it (v.begin ());
170        typename V::iterator it_end (v.end ());
171        while (it != it_end)
172            functor_type::apply (*it, t), ++ it;
173    }
174
175    // Dispatcher
176    template<template <class T1, class T2> class F, class V, class T>
177    BOOST_UBLAS_INLINE
178    void vector_assign_scalar (V &v, const T &t) {
179        typedef typename V::storage_category storage_category;
180        vector_assign_scalar<F> (v, t, storage_category ());
181    }
182
183    template<class SC, bool COMPUTED, class RI>
184    struct vector_assign_traits {
185        typedef SC storage_category;
186    };
187
188    template<bool COMPUTED>
189    struct vector_assign_traits<dense_tag, COMPUTED, packed_random_access_iterator_tag> {
190        typedef packed_tag storage_category;
191    };
192    template<>
193    struct vector_assign_traits<dense_tag, false, sparse_bidirectional_iterator_tag> {
194        typedef sparse_tag storage_category;
195    };
196    template<>
197    struct vector_assign_traits<dense_tag, true, sparse_bidirectional_iterator_tag> {
198        typedef sparse_proxy_tag storage_category;
199    };
200
201    template<bool COMPUTED>
202    struct vector_assign_traits<dense_proxy_tag, COMPUTED, packed_random_access_iterator_tag> {
203        typedef packed_proxy_tag storage_category;
204    };
205    template<>
206    struct vector_assign_traits<dense_proxy_tag, false, sparse_bidirectional_iterator_tag> {
207        typedef sparse_proxy_tag storage_category;
208    };
209    template<>
210    struct vector_assign_traits<dense_proxy_tag, true, sparse_bidirectional_iterator_tag> {
211        typedef sparse_proxy_tag storage_category;
212    };
213
214    template<>
215    struct vector_assign_traits<packed_tag, false, sparse_bidirectional_iterator_tag> {
216        typedef sparse_tag storage_category;
217    };
218    template<>
219    struct vector_assign_traits<packed_tag, true, sparse_bidirectional_iterator_tag> {
220        typedef sparse_proxy_tag storage_category;
221    };
222
223    template<bool COMPUTED>
224    struct vector_assign_traits<packed_proxy_tag, COMPUTED, sparse_bidirectional_iterator_tag> {
225        typedef sparse_proxy_tag storage_category;
226    };
227
228    template<>
229    struct vector_assign_traits<sparse_tag, true, dense_random_access_iterator_tag> {
230        typedef sparse_proxy_tag storage_category;
231    };
232    template<>
233    struct vector_assign_traits<sparse_tag, true, packed_random_access_iterator_tag> {
234        typedef sparse_proxy_tag storage_category;
235    };
236    template<>
237    struct vector_assign_traits<sparse_tag, true, sparse_bidirectional_iterator_tag> {
238        typedef sparse_proxy_tag storage_category;
239    };
240
241    // Explicitly iterating
242    template<template <class T1, class T2> class F, class V, class E>
243    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
244    void iterating_vector_assign (V &v, const vector_expression<E> &e) {
245        typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
246        typedef typename V::difference_type difference_type;
247        difference_type size (BOOST_UBLAS_SAME (v.size (), e ().size ()));
248        typename V::iterator it (v.begin ());
249        BOOST_UBLAS_CHECK (v.end () - it == size, bad_size ());
250        typename E::const_iterator ite (e ().begin ());
251        BOOST_UBLAS_CHECK (e ().end () - ite == size, bad_size ());
252#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
253        while (-- size >= 0)
254            functor_type::apply (*it, *ite), ++ it, ++ ite;
255#else
256        DD (size, 2, r, (functor_type::apply (*it, *ite), ++ it, ++ ite));
257#endif
258    }
259    // Explicitly indexing
260    template<template <class T1, class T2> class F, class V, class E>
261    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
262    void indexing_vector_assign (V &v, const vector_expression<E> &e) {
263        typedef F<typename V::reference, typename E::value_type> functor_type;
264        typedef typename V::size_type size_type;
265        size_type size (BOOST_UBLAS_SAME (v.size (), e ().size ()));
266#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
267        for (size_type i = 0; i < size; ++ i)
268            functor_type::apply (v (i), e () (i));
269#else
270        size_type i (0);
271        DD (size, 2, r, (functor_type::apply (v (i), e () (i)), ++ i));
272#endif
273    }
274
275    // Dense (proxy) case
276    template<template <class T1, class T2> class F, class V, class E>
277    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
278    void vector_assign (V &v, const vector_expression<E> &e, dense_proxy_tag) {
279#ifdef BOOST_UBLAS_USE_INDEXING
280        indexing_vector_assign<F> (v, e);
281#elif BOOST_UBLAS_USE_ITERATING
282        iterating_vector_assign<F> (v, e);
283#else
284        typedef typename V::size_type size_type;
285        size_type size (BOOST_UBLAS_SAME (v.size (), e ().size ()));
286        if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
287            iterating_vector_assign<F> (v, e);
288        else
289            indexing_vector_assign<F> (v, e);
290#endif
291    }
292    // Packed (proxy) case
293    template<template <class T1, class T2> class F, class V, class E>
294    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
295    void vector_assign (V &v, const vector_expression<E> &e, packed_proxy_tag) {
296        BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
297        typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
298        typedef typename V::difference_type difference_type;
299        typedef typename V::value_type value_type;
300#if BOOST_UBLAS_TYPE_CHECK
301        vector<value_type> cv (v.size ());
302        indexing_vector_assign<scalar_assign> (cv, v);
303        indexing_vector_assign<F> (cv, e);
304#endif
305        typename V::iterator it (v.begin ());
306        typename V::iterator it_end (v.end ());
307        typename E::const_iterator ite (e ().begin ());
308        typename E::const_iterator ite_end (e ().end ());
309        difference_type it_size (it_end - it);
310        difference_type ite_size (ite_end - ite);
311        if (it_size > 0 && ite_size > 0) {
312            difference_type size ((std::min) (difference_type (it.index () - ite.index ()), ite_size));
313            if (size > 0) {
314                ite += size;
315                ite_size -= size;
316            }
317        }
318        if (it_size > 0 && ite_size > 0) {
319            difference_type size ((std::min) (difference_type (ite.index () - it.index ()), it_size));
320            if (size > 0) {
321                it_size -= size;
322                if (!functor_type::computed) {
323                    while (-- size >= 0)    // zeroing
324                        functor_type::apply (*it, value_type/*zero*/()), ++ it;
325                } else {
326                    it += size;
327                }
328            }
329        }
330        difference_type size ((std::min) (it_size, ite_size));
331        it_size -= size;
332        ite_size -= size;
333        while (-- size >= 0)
334            functor_type::apply (*it, *ite), ++ it, ++ ite;
335        size = it_size;
336        if (!functor_type::computed) {
337            while (-- size >= 0)    // zeroing
338                functor_type::apply (*it, value_type/*zero*/()), ++ it;
339        } else {
340            it += size;
341        }
342#if BOOST_UBLAS_TYPE_CHECK
343        if (! disable_type_check<bool>::value)
344            BOOST_UBLAS_CHECK (detail::expression_type_check (v, cv), external_logic ());
345#endif
346    }
347    // Sparse case
348    template<template <class T1, class T2> class F, class V, class E>
349    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
350    void vector_assign (V &v, const vector_expression<E> &e, sparse_tag) {
351        BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
352        typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
353        BOOST_STATIC_ASSERT ((!functor_type::computed));
354        typedef typename V::value_type value_type;
355#if BOOST_UBLAS_TYPE_CHECK
356        vector<value_type> cv (v.size ());
357        indexing_vector_assign<scalar_assign> (cv, v);
358        indexing_vector_assign<F> (cv, e);
359#endif
360        v.clear ();
361        typename E::const_iterator ite (e ().begin ());
362        typename E::const_iterator ite_end (e ().end ());
363        while (ite != ite_end) {
364            value_type t (*ite);
365            if (t != value_type/*zero*/())
366                v.insert_element (ite.index (), t);
367            ++ ite;
368        }
369#if BOOST_UBLAS_TYPE_CHECK
370        if (! disable_type_check<bool>::value)
371            BOOST_UBLAS_CHECK (detail::expression_type_check (v, cv), external_logic ());
372#endif
373    }
374    // Sparse proxy or functional case
375    template<template <class T1, class T2> class F, class V, class E>
376    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
377    void vector_assign (V &v, const vector_expression<E> &e, sparse_proxy_tag) {
378        BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
379        typedef F<typename V::iterator::reference, typename E::value_type> functor_type;
380        typedef typename V::size_type size_type;
381        typedef typename V::difference_type difference_type;
382        typedef typename V::value_type value_type;
383        typedef typename V::reference reference;
384#if BOOST_UBLAS_TYPE_CHECK
385        vector<value_type> cv (v.size ());
386        indexing_vector_assign<scalar_assign> (cv, v);
387        indexing_vector_assign<F> (cv, e);
388#endif
389        detail::make_conformant (v, e);
390
391        typename V::iterator it (v.begin ());
392        typename V::iterator it_end (v.end ());
393        typename E::const_iterator ite (e ().begin ());
394        typename E::const_iterator ite_end (e ().end ());
395        if (it != it_end && ite != ite_end) {
396            size_type it_index = it.index (), ite_index = ite.index ();
397            while (true) {
398                difference_type compare = it_index - ite_index;
399                if (compare == 0) {
400                    functor_type::apply (*it, *ite);
401                    ++ it, ++ ite;
402                    if (it != it_end && ite != ite_end) {
403                        it_index = it.index ();
404                        ite_index = ite.index ();
405                    } else
406                        break;
407                } else if (compare < 0) {
408                    if (!functor_type::computed) {
409                        functor_type::apply (*it, value_type/*zero*/());
410                        ++ it;
411                    } else
412                        increment (it, it_end, - compare);
413                    if (it != it_end)
414                        it_index = it.index ();
415                    else
416                        break;
417                } else if (compare > 0) {
418                    increment (ite, ite_end, compare);
419                    if (ite != ite_end)
420                        ite_index = ite.index ();
421                    else
422                        break;
423                }
424            }
425        }
426
427        if (!functor_type::computed) {
428            while (it != it_end) {  // zeroing
429                functor_type::apply (*it, value_type/*zero*/());
430                ++ it;
431            }
432        } else {
433            it = it_end;
434        }
435#if BOOST_UBLAS_TYPE_CHECK
436        if (! disable_type_check<bool>::value)
437            BOOST_UBLAS_CHECK (detail::expression_type_check (v, cv), external_logic ());
438#endif
439    }
440
441    // Dispatcher
442    template<template <class T1, class T2> class F, class V, class E>
443    BOOST_UBLAS_INLINE
444    void vector_assign (V &v, const vector_expression<E> &e) {
445        typedef typename vector_assign_traits<typename V::storage_category,
446                                              F<typename V::reference, typename E::value_type>::computed,
447                                              typename E::const_iterator::iterator_category>::storage_category storage_category;
448        vector_assign<F> (v, e, storage_category ());
449    }
450
451    template<class SC, class RI>
452    struct vector_swap_traits {
453        typedef SC storage_category;
454    };
455
456    template<>
457    struct vector_swap_traits<dense_proxy_tag, sparse_bidirectional_iterator_tag> {
458        typedef sparse_proxy_tag storage_category;
459    };
460
461    template<>
462    struct vector_swap_traits<packed_proxy_tag, sparse_bidirectional_iterator_tag> {
463        typedef sparse_proxy_tag storage_category;
464    };
465
466    // Dense (proxy) case
467    template<template <class T1, class T2> class F, class V, class E>
468    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
469    void vector_swap (V &v, vector_expression<E> &e, dense_proxy_tag) {
470        typedef F<typename V::iterator::reference, typename E::iterator::reference> functor_type;
471        typedef typename V::difference_type difference_type;
472        difference_type size (BOOST_UBLAS_SAME (v.size (), e ().size ()));
473        typename V::iterator it (v.begin ());
474        typename E::iterator ite (e ().begin ());
475        while (-- size >= 0)
476            functor_type::apply (*it, *ite), ++ it, ++ ite;
477    }
478    // Packed (proxy) case
479    template<template <class T1, class T2> class F, class V, class E>
480    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
481    void vector_swap (V &v, vector_expression<E> &e, packed_proxy_tag) {
482        typedef F<typename V::iterator::reference, typename E::iterator::reference> functor_type;
483        typedef typename V::difference_type difference_type;
484        typename V::iterator it (v.begin ());
485        typename V::iterator it_end (v.end ());
486        typename E::iterator ite (e ().begin ());
487        typename E::iterator ite_end (e ().end ());
488        difference_type it_size (it_end - it);
489        difference_type ite_size (ite_end - ite);
490        if (it_size > 0 && ite_size > 0) {
491            difference_type size ((std::min) (difference_type (it.index () - ite.index ()), ite_size));
492            if (size > 0) {
493                ite += size;
494                ite_size -= size;
495            }
496        }
497        if (it_size > 0 && ite_size > 0) {
498            difference_type size ((std::min) (difference_type (ite.index () - it.index ()), it_size));
499            if (size > 0)
500                it_size -= size;
501        }
502        difference_type size ((std::min) (it_size, ite_size));
503        it_size -= size;
504        ite_size -= size;
505        while (-- size >= 0)
506            functor_type::apply (*it, *ite), ++ it, ++ ite;
507    }
508    // Sparse proxy case
509    template<template <class T1, class T2> class F, class V, class E>
510    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
511    void vector_swap (V &v, vector_expression<E> &e, sparse_proxy_tag) {
512        BOOST_UBLAS_CHECK (v.size () == e ().size (), bad_size ());
513        typedef F<typename V::iterator::reference, typename E::iterator::reference> functor_type;
514        typedef typename V::size_type size_type;
515        typedef typename V::difference_type difference_type;
516        typedef typename V::value_type value_type;
517
518        detail::make_conformant (v, e);
519        // FIXME should be a seperate restriction for E
520        detail::make_conformant (e (), v);
521
522        typename V::iterator it (v.begin ());
523        typename V::iterator it_end (v.end ());
524        typename E::iterator ite (e ().begin ());
525        typename E::iterator ite_end (e ().end ());
526        if (it != it_end && ite != ite_end) {
527            size_type it_index = it.index (), ite_index = ite.index ();
528            while (true) {
529                difference_type compare = it_index - ite_index;
530                if (compare == 0) {
531                    functor_type::apply (*it, *ite);
532                    ++ it, ++ ite;
533                    if (it != it_end && ite != ite_end) {
534                        it_index = it.index ();
535                        ite_index = ite.index ();
536                    } else
537                        break;
538                } else if (compare < 0) {
539                    increment (it, it_end, - compare);
540                    if (it != it_end)
541                        it_index = it.index ();
542                    else
543                        break;
544                } else if (compare > 0) {
545                    increment (ite, ite_end, compare);
546                    if (ite != ite_end)
547                        ite_index = ite.index ();
548                    else
549                        break;
550                }
551            }
552        }
553
554#if BOOST_UBLAS_TYPE_CHECK
555        increment (ite, ite_end);
556        increment (it, it_end);
557#endif
558    }
559
560    // Dispatcher
561    template<template <class T1, class T2> class F, class V, class E>
562    BOOST_UBLAS_INLINE
563    void vector_swap (V &v, vector_expression<E> &e) {
564        typedef typename vector_swap_traits<typename V::storage_category,
565                                            typename E::const_iterator::iterator_category>::storage_category storage_category;
566        vector_swap<F> (v, e, storage_category ());
567    }
568
569}}}
570
571#endif
Note: See TracBrowser for help on using the repository browser.