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

Revision 857, 93.7 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_HERMITIAN_H
18#define BOOST_UBLAS_HERMITIAN_H
19
20#include <boost/numeric/ublas/matrix.hpp>
21#include <boost/numeric/ublas/detail/temporary.hpp>
22
23// Iterators based on ideas of Jeremy Siek
24// Hermitian matrices are square. Thanks to Peter Schmitteckert for spotting this.
25
26namespace boost { namespace numeric { namespace ublas {
27
28    template<class M>
29    bool is_hermitian (const M &m) {
30        typedef typename M::size_type size_type;
31
32        if (m.size1 () != m.size2 ())
33            return false;
34        size_type size = BOOST_UBLAS_SAME (m.size1 (), m.size2 ());
35        for (size_type i = 0; i < size; ++ i) {
36            for (size_type j = i; j < size; ++ j) {
37                if (m (i, j) != conj (m (j, i)))
38                    return false;
39            }
40        }
41        return true;
42    }
43
44#ifdef BOOST_UBLAS_STRICT_HERMITIAN
45
46    template<class M>
47    class hermitian_matrix_element:
48       public container_reference<M> {
49    public:
50        typedef M matrix_type;
51        typedef typename M::size_type size_type;
52        typedef typename M::value_type value_type;
53        typedef const value_type &const_reference;
54        typedef value_type &reference;
55        typedef value_type *pointer;
56
57        // Construction and destruction
58        BOOST_UBLAS_INLINE
59        hermitian_matrix_element (matrix_type &m, size_type i, size_type j, value_type d):
60            container_reference<matrix_type> (m), i_ (i), j_ (j), d_ (d), dirty_ (false) {}
61        BOOST_UBLAS_INLINE
62        hermitian_matrix_element (const hermitian_matrix_element &p):
63            container_reference<matrix_type> (p), i_ (p.i_), d_ (p.d_), dirty_ (p.dirty_) {}
64        BOOST_UBLAS_INLINE
65        ~hermitian_matrix_element () {
66            if (dirty_)
67                (*this) ().insert_element (i_, j_, d_);
68        }
69
70        // Assignment
71        BOOST_UBLAS_INLINE
72        hermitian_matrix_element &operator = (const hermitian_matrix_element &p) {
73            // Overide the implict copy assignment
74            d_ = p.d_;
75            dirty_ = true;
76            return *this;
77        }
78        template<class D>
79        BOOST_UBLAS_INLINE
80        hermitian_matrix_element &operator = (const D &d) {
81            d_ = d;
82            dirty_ = true;
83            return *this;
84        }
85        template<class D>
86        BOOST_UBLAS_INLINE
87        hermitian_matrix_element &operator += (const D &d) {
88            d_ += d;
89            dirty_ = true;
90            return *this;
91        }
92        template<class D>
93        BOOST_UBLAS_INLINE
94        hermitian_matrix_element &operator -= (const D &d) {
95            d_ -= d;
96            dirty_ = true;
97            return *this;
98        }
99        template<class D>
100        BOOST_UBLAS_INLINE
101        hermitian_matrix_element &operator *= (const D &d) {
102            d_ *= d;
103            dirty_ = true;
104            return *this;
105        }
106        template<class D>
107        BOOST_UBLAS_INLINE
108        hermitian_matrix_element &operator /= (const D &d) {
109            d_ /= d;
110            dirty_ = true;
111            return *this;
112        }
113       
114        // Comparison
115        template<class D>
116        BOOST_UBLAS_INLINE
117        bool operator == (const D &d) const {
118            return d_ == d;
119        }
120        template<class D>
121        BOOST_UBLAS_INLINE
122        bool operator != (const D &d) const {
123            return d_ != d;
124        }
125
126        // Conversion
127        BOOST_UBLAS_INLINE
128        operator const_reference () const {
129            return d_;
130        }
131
132        // Swapping
133        BOOST_UBLAS_INLINE
134        void swap (hermitian_matrix_element p) {
135            if (this != &p) {
136                dirty_ = true;
137                p.dirty_ = true;
138                std::swap (d_, p.d_);
139            }
140        }
141        BOOST_UBLAS_INLINE
142        friend void swap (hermitian_matrix_element p1, hermitian_matrix_element p2) {
143            p1.swap (p2);
144        }
145
146    private:
147        size_type i_;
148        size_type j_;
149        value_type d_;
150        bool dirty_;
151    };
152
153    template<class M>
154    struct type_traits<hermitian_matrix_element<M> > {
155        typedef typename M::value_type element_type;
156        typedef type_traits<hermitian_matrix_element<M> > self_type;
157        typedef typename type_traits<element_type>::value_type value_type;
158        typedef typename type_traits<element_type>::const_reference const_reference;
159        typedef hermitian_matrix_element<M> reference;
160        typedef typename type_traits<element_type>::real_type real_type;
161        typedef typename type_traits<element_type>::precision_type precision_type;
162
163        static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
164        static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
165
166        static
167        BOOST_UBLAS_INLINE
168        real_type real (const_reference t) {
169            return type_traits<element_type>::real (t);
170        }
171        static
172        BOOST_UBLAS_INLINE
173        real_type imag (const_reference t) {
174            return type_traits<element_type>::imag (t);
175        }
176        static
177        BOOST_UBLAS_INLINE
178        value_type conj (const_reference t) {
179            return type_traits<element_type>::conj (t);
180        }
181
182        static
183        BOOST_UBLAS_INLINE
184        real_type abs (const_reference t) {
185            return type_traits<element_type>::abs (t);
186        }
187        static
188        BOOST_UBLAS_INLINE
189        value_type sqrt (const_reference t) {
190            return type_traits<element_type>::sqrt (t);
191        }
192
193        static
194        BOOST_UBLAS_INLINE
195        real_type norm_1 (const_reference t) {
196            return type_traits<element_type>::norm_1 (t);
197        }
198        static
199        BOOST_UBLAS_INLINE
200        real_type norm_2 (const_reference t) {
201            return type_traits<element_type>::norm_2 (t);
202        }
203        static
204        BOOST_UBLAS_INLINE
205        real_type norm_inf (const_reference t) {
206            return type_traits<element_type>::norm_inf (t);
207        }
208
209        static
210        BOOST_UBLAS_INLINE
211        bool equals (const_reference t1, const_reference t2) {
212            return type_traits<element_type>::equals (t1, t2);
213        }
214    };
215
216    template<class M1, class T2>
217    struct promote_traits<hermitian_matrix_element<M1>, T2> {
218        typedef typename promote_traits<typename hermitian_matrix_element<M1>::value_type, T2>::promote_type promote_type;
219    };
220    template<class T1, class M2>
221    struct promote_traits<T1, hermitian_matrix_element<M2> > {
222        typedef typename promote_traits<T1, typename hermitian_matrix_element<M2>::value_type>::promote_type promote_type;
223    };
224    template<class M1, class M2>
225    struct promote_traits<hermitian_matrix_element<M1>, hermitian_matrix_element<M2> > {
226        typedef typename promote_traits<typename hermitian_matrix_element<M1>::value_type,
227                                        typename hermitian_matrix_element<M2>::value_type>::promote_type promote_type;
228    };
229
230#endif
231
232    // Array based hermitian matrix class
233    template<class T, class TRI, class L, class A>
234    class hermitian_matrix:
235        public matrix_container<hermitian_matrix<T, TRI, L, A> > {
236
237        typedef T &true_reference;
238        typedef T *pointer;
239        typedef TRI triangular_type;
240        typedef L layout_type;
241        typedef hermitian_matrix<T, TRI, L, A> self_type;
242    public:
243#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
244        using matrix_container<self_type>::operator ();
245#endif
246        typedef typename A::size_type size_type;
247        typedef typename A::difference_type difference_type;
248        typedef T value_type;
249        // FIXME no better way to not return the address of a temporary?
250        // typedef const T &const_reference;
251        typedef const T const_reference;
252#ifndef BOOST_UBLAS_STRICT_HERMITIAN
253        typedef T &reference;
254#else
255        typedef hermitian_matrix_element<self_type> reference;
256#endif
257        typedef A array_type;
258
259        typedef const matrix_reference<const self_type> const_closure_type;
260        typedef matrix_reference<self_type> closure_type;
261        typedef vector<T, A> vector_temporary_type;
262        typedef matrix<T, L, A> matrix_temporary_type;  // general sub-matrix
263        typedef packed_tag storage_category;
264        typedef typename L::orientation_category orientation_category;
265
266        // Construction and destruction
267        BOOST_UBLAS_INLINE
268        hermitian_matrix ():
269            matrix_container<self_type> (),
270            size_ (0), data_ (0) {}
271        BOOST_UBLAS_INLINE
272        hermitian_matrix (size_type size):
273            matrix_container<self_type> (),
274            size_ (BOOST_UBLAS_SAME (size, size)), data_ (triangular_type::packed_size (layout_type (), size, size)) {
275        }
276        BOOST_UBLAS_INLINE
277        hermitian_matrix (size_type size1, size_type size2):
278            matrix_container<self_type> (),
279            size_ (BOOST_UBLAS_SAME (size1, size2)), data_ (triangular_type::packed_size (layout_type (), size1, size2)) {
280        }
281        BOOST_UBLAS_INLINE
282        hermitian_matrix (size_type size, const array_type &data):
283            matrix_container<self_type> (),
284            size_ (size), data_ (data) {}
285        BOOST_UBLAS_INLINE
286        hermitian_matrix (const hermitian_matrix &m):
287            matrix_container<self_type> (),
288            size_ (m.size_), data_ (m.data_) {}
289        template<class AE>
290        BOOST_UBLAS_INLINE
291        hermitian_matrix (const matrix_expression<AE> &ae):
292            matrix_container<self_type> (),
293            size_ (BOOST_UBLAS_SAME (ae ().size1 (), ae ().size2 ())),
294            data_ (triangular_type::packed_size (layout_type (), size_, size_)) {
295            matrix_assign<scalar_assign> (*this, ae);
296        }
297
298        // Accessors
299        BOOST_UBLAS_INLINE
300        size_type size1 () const {
301            return size_;
302        }
303        BOOST_UBLAS_INLINE
304        size_type size2 () const {
305            return size_;
306        }
307
308        // Storage accessors
309        BOOST_UBLAS_INLINE
310        const array_type &data () const {
311            return data_;
312        }
313        BOOST_UBLAS_INLINE
314        array_type &data () {
315            return data_;
316        }
317
318        // Resizing
319        BOOST_UBLAS_INLINE
320        void resize (size_type size, bool preserve = true) {
321            if (preserve) {
322                self_type temporary (size, size);
323                detail::matrix_resize_preserve<layout_type> (*this, temporary);
324            }
325            else {
326                data ().resize (triangular_type::packed_size (layout_type (), size, size));
327                size_ = size;
328            }
329        }
330        BOOST_UBLAS_INLINE
331        void resize (size_type size1, size_type size2, bool preserve = true) {
332            resize (BOOST_UBLAS_SAME (size1, size2), preserve);
333        }
334        BOOST_UBLAS_INLINE
335        void resize_packed_preserve (size_type size) {
336            size_ = BOOST_UBLAS_SAME (size, size);
337            data ().resize (triangular_type::packed_size (layout_type (), size_, size_), value_type ());
338        }
339
340        // Element access
341        BOOST_UBLAS_INLINE
342        const_reference operator () (size_type i, size_type j) const {
343            BOOST_UBLAS_CHECK (i < size_, bad_index ());
344            BOOST_UBLAS_CHECK (j < size_, bad_index ());
345            // if (i == j)
346            //    return type_traits<value_type>::real (data () [triangular_type::element (layout_type (), i, size_, i, size_)]);
347            // else
348            if (triangular_type::other (i, j))
349                return data () [triangular_type::element (layout_type (), i, size_, j, size_)];
350            else
351                return type_traits<value_type>::conj (data () [triangular_type::element (layout_type (), j, size_, i, size_)]);
352        }
353        BOOST_UBLAS_INLINE
354        true_reference at_element (size_type i, size_type j) {
355            BOOST_UBLAS_CHECK (i < size_, bad_index ());
356            BOOST_UBLAS_CHECK (j < size_, bad_index ());
357            BOOST_UBLAS_CHECK (triangular_type::other (i, j), bad_index ());
358            return data () [triangular_type::element (layout_type (), i, size_, j, size_)];
359        }
360        BOOST_UBLAS_INLINE
361        reference operator () (size_type i, size_type j) {
362#ifndef BOOST_UBLAS_STRICT_HERMITIAN
363            if (triangular_type::other (i, j))
364                return at_element (i, j);
365            else {
366                external_logic ().raise ();
367                // arbitary return value
368                return data () [triangular_type::element (layout_type (), j, size_, i, size_)];
369            }
370#else
371        if (triangular_type::other (i, j))
372            return reference (*this, i, j, data () [triangular_type::element (layout_type (), i, size_, j, size_)]);
373        else
374            return reference (*this, i, j, type_traits<value_type>::conj (data () [triangular_type::element (layout_type (), j, size_, i, size_)]));
375#endif
376        }
377
378        // Element assignemnt
379        BOOST_UBLAS_INLINE
380        true_reference insert_element (size_type i, size_type j, const_reference t) {
381            BOOST_UBLAS_CHECK (i < size_, bad_index ());
382            BOOST_UBLAS_CHECK (j < size_, bad_index ());
383            if (triangular_type::other (i, j)) {
384                return (data () [triangular_type::element (layout_type (), i, size_, j, size_)] = t);
385            } else {
386                return (data () [triangular_type::element (layout_type (), j, size_, i, size_)] = type_traits<value_type>::conj (t));
387            }
388        }
389        BOOST_UBLAS_INLINE
390        void erase_element (size_type i, size_type j) {
391            BOOST_UBLAS_CHECK (i < size_, bad_index ());
392            BOOST_UBLAS_CHECK (j < size_, bad_index ());
393            data () [triangular_type::element (layout_type (), i, size_, j, size_)] = value_type/*zero*/();
394        }
395
396        // Zeroing
397        BOOST_UBLAS_INLINE
398        void clear () {
399            std::fill (data ().begin (), data ().end (), value_type/*zero*/());
400        }
401
402        // Assignment
403        BOOST_UBLAS_INLINE
404        hermitian_matrix &operator = (const hermitian_matrix &m) {
405            size_ = m.size_;
406            data () = m.data ();
407            return *this;
408        }
409        BOOST_UBLAS_INLINE
410        hermitian_matrix &assign_temporary (hermitian_matrix &m) {
411            swap (m);
412            return *this;
413        }
414        template<class AE>
415        BOOST_UBLAS_INLINE
416        hermitian_matrix &operator = (const matrix_expression<AE> &ae) {
417            self_type temporary (ae);
418            return assign_temporary (temporary);
419        }
420        template<class AE>
421        BOOST_UBLAS_INLINE
422        hermitian_matrix &assign (const matrix_expression<AE> &ae) {
423            matrix_assign<scalar_assign> (*this, ae);
424            return *this;
425        }
426        template<class AE>
427        BOOST_UBLAS_INLINE
428        hermitian_matrix& operator += (const matrix_expression<AE> &ae) {
429            self_type temporary (*this + ae);
430            return assign_temporary (temporary);
431        }
432        template<class AE>
433        BOOST_UBLAS_INLINE
434        hermitian_matrix &plus_assign (const matrix_expression<AE> &ae) {
435            matrix_assign<scalar_plus_assign> (*this, ae);
436            return *this;
437        }
438        template<class AE>
439        BOOST_UBLAS_INLINE
440        hermitian_matrix& operator -= (const matrix_expression<AE> &ae) {
441            self_type temporary (*this - ae);
442            return assign_temporary (temporary);
443        }
444        template<class AE>
445        BOOST_UBLAS_INLINE
446        hermitian_matrix &minus_assign (const matrix_expression<AE> &ae) {
447            matrix_assign<scalar_minus_assign> (*this, ae);
448            return *this;
449        }
450        template<class AT>
451        BOOST_UBLAS_INLINE
452        hermitian_matrix& operator *= (const AT &at) {
453            // Multiplication is only allowed for real scalars,
454            // otherwise the resulting matrix isn't hermitian.
455            // Thanks to Peter Schmitteckert for spotting this.
456            BOOST_UBLAS_CHECK (type_traits<value_type>::imag (at) == 0, non_real ());
457            matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
458            return *this;
459        }
460        template<class AT>
461        BOOST_UBLAS_INLINE
462        hermitian_matrix& operator /= (const AT &at) {
463            // Multiplication is only allowed for real scalars,
464            // otherwise the resulting matrix isn't hermitian.
465            // Thanks to Peter Schmitteckert for spotting this.
466            BOOST_UBLAS_CHECK (type_traits<value_type>::imag (at) == 0, non_real ());
467            matrix_assign_scalar<scalar_divides_assign> (*this, at);
468            return *this;
469        }
470
471        // Swapping
472        BOOST_UBLAS_INLINE
473        void swap (hermitian_matrix &m) {
474            if (this != &m) {
475                std::swap (size_, m.size_);
476                data ().swap (m.data ());
477            }
478        }
479        BOOST_UBLAS_INLINE
480        friend void swap (hermitian_matrix &m1, hermitian_matrix &m2) {
481            m1.swap (m2);
482        }
483
484        // Iterator types
485#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
486        typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
487        typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
488        typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
489        typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
490#else
491        class const_iterator1;
492        class iterator1;
493        class const_iterator2;
494        class iterator2;
495#endif
496        typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
497        typedef reverse_iterator_base1<iterator1> reverse_iterator1;
498        typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
499        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
500
501        // Element lookup
502        BOOST_UBLAS_INLINE
503        const_iterator1 find1 (int /* rank */, size_type i, size_type j) const {
504            return const_iterator1 (*this, i, j);
505        }
506        BOOST_UBLAS_INLINE
507        iterator1 find1 (int rank, size_type i, size_type j) {
508            if (rank == 1)
509                i = triangular_type::mutable_restrict1 (i, j);
510            return iterator1 (*this, i, j);
511        }
512        BOOST_UBLAS_INLINE
513        const_iterator2 find2 (int /* rank */, size_type i, size_type j) const {
514            return const_iterator2 (*this, i, j);
515        }
516        BOOST_UBLAS_INLINE
517        iterator2 find2 (int rank, size_type i, size_type j) {
518            if (rank == 1)
519                j = triangular_type::mutable_restrict2 (i, j);
520            return iterator2 (*this, i, j);
521        }
522
523        // Iterators simply are indices.
524
525#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
526        class const_iterator1:
527            public container_const_reference<hermitian_matrix>,
528            public random_access_iterator_base<packed_random_access_iterator_tag,
529                                               const_iterator1, value_type> {
530        public:
531            typedef typename hermitian_matrix::value_type value_type;
532            typedef typename hermitian_matrix::difference_type difference_type;
533            typedef typename hermitian_matrix::const_reference reference;
534            typedef const typename hermitian_matrix::pointer pointer;
535
536            typedef const_iterator2 dual_iterator_type;
537            typedef const_reverse_iterator2 dual_reverse_iterator_type;
538
539            // Construction and destruction
540            BOOST_UBLAS_INLINE
541            const_iterator1 ():
542                container_const_reference<self_type> (), it1_ (), it2_ () {}
543            BOOST_UBLAS_INLINE
544            const_iterator1 (const self_type &m, size_type it1, size_type it2):
545                container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
546            BOOST_UBLAS_INLINE
547            const_iterator1 (const iterator1 &it):
548                container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
549
550            // Arithmetic
551            BOOST_UBLAS_INLINE
552            const_iterator1 &operator ++ () {
553                ++ it1_;
554                return *this;
555            }
556            BOOST_UBLAS_INLINE
557            const_iterator1 &operator -- () {
558                -- it1_;
559                return *this;
560            }
561            BOOST_UBLAS_INLINE
562            const_iterator1 &operator += (difference_type n) {
563                it1_ += n;
564                return *this;
565            }
566            BOOST_UBLAS_INLINE
567            const_iterator1 &operator -= (difference_type n) {
568                it1_ -= n;
569                return *this;
570            }
571            BOOST_UBLAS_INLINE
572            difference_type operator - (const const_iterator1 &it) const {
573                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
574                BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
575                return it1_ - it.it1_;
576            }
577
578            // Dereference
579            BOOST_UBLAS_INLINE
580            const_reference operator * () const {
581                return (*this) () (it1_, it2_);
582            }
583
584#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
585            BOOST_UBLAS_INLINE
586#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
587            typename self_type::
588#endif
589            const_iterator2 begin () const {
590                return (*this) ().find2 (1, it1_, 0);
591            }
592            BOOST_UBLAS_INLINE
593#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
594            typename self_type::
595#endif
596            const_iterator2 end () const {
597                return (*this) ().find2 (1, it1_, (*this) ().size2 ());
598            }
599            BOOST_UBLAS_INLINE
600#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
601            typename self_type::
602#endif
603            const_reverse_iterator2 rbegin () const {
604                return const_reverse_iterator2 (end ());
605            }
606            BOOST_UBLAS_INLINE
607#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
608            typename self_type::
609#endif
610            const_reverse_iterator2 rend () const {
611                return const_reverse_iterator2 (begin ());
612            }
613#endif
614
615            // Indices
616            BOOST_UBLAS_INLINE
617            size_type index1 () const {
618                return it1_;
619            }
620            BOOST_UBLAS_INLINE
621            size_type index2 () const {
622                return it2_;
623            }
624
625            // Assignment
626            BOOST_UBLAS_INLINE
627            const_iterator1 &operator = (const const_iterator1 &it) {
628                container_const_reference<self_type>::assign (&it ());
629                it1_ = it.it1_;
630                it2_ = it.it2_;
631                return *this;
632            }
633
634            // Comparison
635            BOOST_UBLAS_INLINE
636            bool operator == (const const_iterator1 &it) const {
637                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
638                BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
639                return it1_ == it.it1_;
640            }
641            BOOST_UBLAS_INLINE
642            bool operator < (const const_iterator1 &it) const {
643                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
644                BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
645                return it1_ < it.it1_;
646            }
647
648        private:
649            size_type it1_;
650            size_type it2_;
651        };
652#endif
653
654        BOOST_UBLAS_INLINE
655        const_iterator1 begin1 () const {
656            return find1 (0, 0, 0);
657        }
658        BOOST_UBLAS_INLINE
659        const_iterator1 end1 () const {
660            return find1 (0, size_, 0);
661        }
662
663#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
664        class iterator1:
665            public container_reference<hermitian_matrix>,
666            public random_access_iterator_base<packed_random_access_iterator_tag,
667                                               iterator1, value_type> {
668        public:
669            typedef typename hermitian_matrix::value_type value_type;
670            typedef typename hermitian_matrix::difference_type difference_type;
671            typedef typename hermitian_matrix::true_reference reference;
672            typedef typename hermitian_matrix::pointer pointer;
673
674            typedef iterator2 dual_iterator_type;
675            typedef reverse_iterator2 dual_reverse_iterator_type;
676
677            // Construction and destruction
678            BOOST_UBLAS_INLINE
679            iterator1 ():
680                container_reference<self_type> (), it1_ (), it2_ () {}
681            BOOST_UBLAS_INLINE
682            iterator1 (self_type &m, size_type it1, size_type it2):
683                container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
684
685            // Arithmetic
686            BOOST_UBLAS_INLINE
687            iterator1 &operator ++ () {
688                ++ it1_;
689                return *this;
690            }
691            BOOST_UBLAS_INLINE
692            iterator1 &operator -- () {
693                -- it1_;
694                return *this;
695            }
696            BOOST_UBLAS_INLINE
697            iterator1 &operator += (difference_type n) {
698                it1_ += n;
699                return *this;
700            }
701            BOOST_UBLAS_INLINE
702            iterator1 &operator -= (difference_type n) {
703                it1_ -= n;
704                return *this;
705            }
706            BOOST_UBLAS_INLINE
707            difference_type operator - (const iterator1 &it) const {
708                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
709                BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
710                return it1_ - it.it1_;
711            }
712
713            // Dereference
714            BOOST_UBLAS_INLINE
715            reference operator * () const {
716                return (*this) ().at_element (it1_, it2_);
717            }
718
719#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
720            BOOST_UBLAS_INLINE
721#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
722            typename self_type::
723#endif
724            iterator2 begin () const {
725                return (*this) ().find2 (1, it1_, 0);
726            }
727            BOOST_UBLAS_INLINE
728#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
729            typename self_type::
730#endif
731            iterator2 end () const {
732                return (*this) ().find2 (1, it1_, (*this) ().size2 ());
733            }
734            BOOST_UBLAS_INLINE
735#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
736            typename self_type::
737#endif
738            reverse_iterator2 rbegin () const {
739                return reverse_iterator2 (end ());
740            }
741            BOOST_UBLAS_INLINE
742#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
743            typename self_type::
744#endif
745            reverse_iterator2 rend () const {
746                return reverse_iterator2 (begin ());
747            }
748#endif
749
750            // Indices
751            BOOST_UBLAS_INLINE
752            size_type index1 () const {
753                return it1_;
754            }
755            BOOST_UBLAS_INLINE
756            size_type index2 () const {
757                return it2_;
758            }
759
760            // Assignment
761            BOOST_UBLAS_INLINE
762            iterator1 &operator = (const iterator1 &it) {
763                container_reference<self_type>::assign (&it ());
764                it1_ = it.it1_;
765                it2_ = it.it2_;
766                return *this;
767            }
768
769            // Comparison
770            BOOST_UBLAS_INLINE
771            bool operator == (const iterator1 &it) const {
772                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
773                BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
774                return it1_ == it.it1_;
775            }
776            BOOST_UBLAS_INLINE
777            bool operator < (const iterator1 &it) const {
778                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
779                BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
780                return it1_ < it.it1_;
781            }
782
783        private:
784            size_type it1_;
785            size_type it2_;
786
787            friend class const_iterator1;
788        };
789#endif
790
791        BOOST_UBLAS_INLINE
792        iterator1 begin1 () {
793            return find1 (0, 0, 0);
794        }
795        BOOST_UBLAS_INLINE
796        iterator1 end1 () {
797            return find1 (0, size_, 0);
798        }
799
800#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
801        class const_iterator2:
802            public container_const_reference<hermitian_matrix>,
803            public random_access_iterator_base<packed_random_access_iterator_tag,
804                                               const_iterator2, value_type> {
805        public:
806            typedef typename hermitian_matrix::value_type value_type;
807            typedef typename hermitian_matrix::difference_type difference_type;
808            typedef typename hermitian_matrix::const_reference reference;
809            typedef const typename hermitian_matrix::pointer pointer;
810
811            typedef const_iterator1 dual_iterator_type;
812            typedef const_reverse_iterator1 dual_reverse_iterator_type;
813
814            // Construction and destruction
815            BOOST_UBLAS_INLINE
816            const_iterator2 ():
817                container_const_reference<self_type> (), it1_ (), it2_ () {}
818            BOOST_UBLAS_INLINE
819            const_iterator2 (const self_type &m, size_type it1, size_type it2):
820                container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
821            BOOST_UBLAS_INLINE
822            const_iterator2 (const iterator2 &it):
823                container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
824
825            // Arithmetic
826            BOOST_UBLAS_INLINE
827            const_iterator2 &operator ++ () {
828                ++ it2_;
829                return *this;
830            }
831            BOOST_UBLAS_INLINE
832            const_iterator2 &operator -- () {
833                -- it2_;
834                return *this;
835            }
836            BOOST_UBLAS_INLINE
837            const_iterator2 &operator += (difference_type n) {
838                it2_ += n;
839                return *this;
840            }
841            BOOST_UBLAS_INLINE
842            const_iterator2 &operator -= (difference_type n) {
843                it2_ -= n;
844                return *this;
845            }
846            BOOST_UBLAS_INLINE
847            difference_type operator - (const const_iterator2 &it) const {
848                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
849                BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
850                return it2_ - it.it2_;
851            }
852
853            // Dereference
854            BOOST_UBLAS_INLINE
855            const_reference operator * () const {
856                return (*this) () (it1_, it2_);
857            }
858
859#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
860            BOOST_UBLAS_INLINE
861#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
862            typename self_type::
863#endif
864            const_iterator1 begin () const {
865                return (*this) ().find1 (1, 0, it2_);
866            }
867            BOOST_UBLAS_INLINE
868#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
869            typename self_type::
870#endif
871            const_iterator1 end () const {
872                return (*this) ().find1 (1, (*this) ().size1 (), it2_);
873            }
874            BOOST_UBLAS_INLINE
875#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
876            typename self_type::
877#endif
878            const_reverse_iterator1 rbegin () const {
879                return const_reverse_iterator1 (end ());
880            }
881            BOOST_UBLAS_INLINE
882#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
883            typename self_type::
884#endif
885            const_reverse_iterator1 rend () const {
886                return const_reverse_iterator1 (begin ());
887            }
888#endif
889
890            // Indices
891            BOOST_UBLAS_INLINE
892            size_type index1 () const {
893                return it1_;
894            }
895            BOOST_UBLAS_INLINE
896            size_type index2 () const {
897                return it2_;
898            }
899
900            // Assignment
901            BOOST_UBLAS_INLINE
902            const_iterator2 &operator = (const const_iterator2 &it) {
903                container_const_reference<self_type>::assign (&it ());
904                it1_ = it.it1_;
905                it2_ = it.it2_;
906                return *this;
907            }
908
909            // Comparison
910            BOOST_UBLAS_INLINE
911            bool operator == (const const_iterator2 &it) const {
912                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
913                BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
914                return it2_ == it.it2_;
915            }
916            BOOST_UBLAS_INLINE
917            bool operator < (const const_iterator2 &it) const {
918                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
919                BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
920                return it2_ < it.it2_;
921            }
922
923        private:
924            size_type it1_;
925            size_type it2_;
926        };
927#endif
928
929        BOOST_UBLAS_INLINE
930        const_iterator2 begin2 () const {
931            return find2 (0, 0, 0);
932        }
933        BOOST_UBLAS_INLINE
934        const_iterator2 end2 () const {
935            return find2 (0, 0, size_);
936        }
937
938#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
939        class iterator2:
940            public container_reference<hermitian_matrix>,
941            public random_access_iterator_base<packed_random_access_iterator_tag,
942                                               iterator2, value_type> {
943        public:
944            typedef typename hermitian_matrix::value_type value_type;
945            typedef typename hermitian_matrix::difference_type difference_type;
946            typedef typename hermitian_matrix::true_reference reference;
947            typedef typename hermitian_matrix::pointer pointer;
948
949            typedef iterator1 dual_iterator_type;
950            typedef reverse_iterator1 dual_reverse_iterator_type;
951
952            // Construction and destruction
953            BOOST_UBLAS_INLINE
954            iterator2 ():
955                container_reference<self_type> (), it1_ (), it2_ () {}
956            BOOST_UBLAS_INLINE
957            iterator2 (self_type &m, size_type it1, size_type it2):
958                container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
959
960            // Arithmetic
961            BOOST_UBLAS_INLINE
962            iterator2 &operator ++ () {
963                ++ it2_;
964                return *this;
965            }
966            BOOST_UBLAS_INLINE
967            iterator2 &operator -- () {
968                -- it2_;
969                return *this;
970            }
971            BOOST_UBLAS_INLINE
972            iterator2 &operator += (difference_type n) {
973                it2_ += n;
974                return *this;
975            }
976            BOOST_UBLAS_INLINE
977            iterator2 &operator -= (difference_type n) {
978                it2_ -= n;
979                return *this;
980            }
981            BOOST_UBLAS_INLINE
982            difference_type operator - (const iterator2 &it) const {
983                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
984                BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
985                return it2_ - it.it2_;
986            }
987
988            // Dereference
989            BOOST_UBLAS_INLINE
990            reference operator * () const {
991                return (*this) ().at_element (it1_, it2_);
992            }
993
994#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
995            BOOST_UBLAS_INLINE
996#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
997            typename self_type::
998#endif
999            iterator1 begin () const {
1000                return (*this) ().find1 (1, 0, it2_);
1001            }
1002            BOOST_UBLAS_INLINE
1003#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1004            typename self_type::
1005#endif
1006            iterator1 end () const {
1007                return (*this) ().find1 (1, (*this) ().size1 (), it2_);
1008            }
1009            BOOST_UBLAS_INLINE
1010#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1011            typename self_type::
1012#endif
1013            reverse_iterator1 rbegin () const {
1014                return reverse_iterator1 (end ());
1015            }
1016            BOOST_UBLAS_INLINE
1017#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1018            typename self_type::
1019#endif
1020            reverse_iterator1 rend () const {
1021                return reverse_iterator1 (begin ());
1022            }
1023#endif
1024
1025            // Indices
1026            BOOST_UBLAS_INLINE
1027            size_type index1 () const {
1028                return it1_;
1029            }
1030            BOOST_UBLAS_INLINE
1031            size_type index2 () const {
1032                return it2_;
1033            }
1034
1035            // Assignment
1036            BOOST_UBLAS_INLINE
1037            iterator2 &operator = (const iterator2 &it) {
1038                container_reference<self_type>::assign (&it ());
1039                it1_ = it.it1_;
1040                it2_ = it.it2_;
1041                return *this;
1042            }
1043
1044            // Comparison
1045            BOOST_UBLAS_INLINE
1046            bool operator == (const iterator2 &it) const {
1047                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1048                BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1049                return it2_ == it.it2_;
1050            }
1051            BOOST_UBLAS_INLINE
1052            bool operator < (const iterator2 &it) const {
1053                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1054                BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1055                return it2_ < it.it2_;
1056            }
1057
1058        private:
1059            size_type it1_;
1060            size_type it2_;
1061
1062            friend class const_iterator2;
1063        };
1064#endif
1065
1066        BOOST_UBLAS_INLINE
1067        iterator2 begin2 () {
1068            return find2 (0, 0, 0);
1069        }
1070        BOOST_UBLAS_INLINE
1071        iterator2 end2 () {
1072            return find2 (0, 0, size_);
1073        }
1074
1075        // Reverse iterators
1076
1077        BOOST_UBLAS_INLINE
1078        const_reverse_iterator1 rbegin1 () const {
1079            return const_reverse_iterator1 (end1 ());
1080        }
1081        BOOST_UBLAS_INLINE
1082        const_reverse_iterator1 rend1 () const {
1083            return const_reverse_iterator1 (begin1 ());
1084        }
1085
1086        BOOST_UBLAS_INLINE
1087        reverse_iterator1 rbegin1 () {
1088            return reverse_iterator1 (end1 ());
1089        }
1090        BOOST_UBLAS_INLINE
1091        reverse_iterator1 rend1 () {
1092            return reverse_iterator1 (begin1 ());
1093        }
1094
1095        BOOST_UBLAS_INLINE
1096        const_reverse_iterator2 rbegin2 () const {
1097            return const_reverse_iterator2 (end2 ());
1098        }
1099        BOOST_UBLAS_INLINE
1100        const_reverse_iterator2 rend2 () const {
1101            return const_reverse_iterator2 (begin2 ());
1102        }
1103
1104        BOOST_UBLAS_INLINE
1105        reverse_iterator2 rbegin2 () {
1106            return reverse_iterator2 (end2 ());
1107        }
1108        BOOST_UBLAS_INLINE
1109        reverse_iterator2 rend2 () {
1110            return reverse_iterator2 (begin2 ());
1111        }
1112
1113    private:
1114        size_type size_;
1115        array_type data_;
1116    };
1117
1118
1119    // Hermitian matrix adaptor class
1120    template<class M, class TRI>
1121    class hermitian_adaptor:
1122        public matrix_expression<hermitian_adaptor<M, TRI> > {
1123
1124        typedef hermitian_adaptor<M, TRI> self_type;
1125        typedef typename M::value_type &true_reference;
1126    public:
1127#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1128        using matrix_expression<self_type>::operator ();
1129#endif
1130        typedef const M const_matrix_type;
1131        typedef M matrix_type;
1132        typedef TRI triangular_type;
1133        typedef typename M::size_type size_type;
1134        typedef typename M::difference_type difference_type;
1135        typedef typename M::value_type value_type;
1136        typedef typename M::value_type const_reference;
1137#ifndef BOOST_UBLAS_STRICT_HERMITIAN
1138        typedef typename boost::mpl::if_<boost::is_const<M>,
1139                                          typename M::value_type,
1140                                          typename M::reference>::type reference;
1141#else
1142        typedef typename boost::mpl::if_<boost::is_const<M>,
1143                                          typename M::value_type,
1144                                          hermitian_matrix_element<self_type> >::type reference;
1145#endif
1146        typedef typename boost::mpl::if_<boost::is_const<M>,
1147                                          typename M::const_closure_type,
1148                                          typename M::closure_type>::type matrix_closure_type;
1149        typedef const self_type const_closure_type;
1150        typedef self_type closure_type;
1151        // Replaced by _temporary_traits to avoid type requirements on M
1152        //typedef typename M::vector_temporary_type vector_temporary_type;
1153        //typedef typename M::matrix_temporary_type matrix_temporary_type;
1154        typedef typename storage_restrict_traits<typename M::storage_category,
1155                                                 packed_proxy_tag>::storage_category storage_category;
1156        typedef typename M::orientation_category orientation_category;
1157
1158        // Construction and destruction
1159        BOOST_UBLAS_INLINE
1160        hermitian_adaptor (matrix_type &data):
1161            matrix_expression<self_type> (),
1162            data_ (data) {
1163            BOOST_UBLAS_CHECK (data_.size1 () == data_.size2 (), bad_size ());
1164        }
1165        BOOST_UBLAS_INLINE
1166        hermitian_adaptor (const hermitian_adaptor &m):
1167            matrix_expression<self_type> (),
1168            data_ (m.data_) {
1169            BOOST_UBLAS_CHECK (data_.size1 () == data_.size2 (), bad_size ());
1170        }
1171
1172        // Accessors
1173        BOOST_UBLAS_INLINE
1174        size_type size1 () const {
1175            return data_.size1 ();
1176        }
1177        BOOST_UBLAS_INLINE
1178        size_type size2 () const {
1179            return data_.size2 ();
1180        }
1181
1182        // Storage accessors
1183        BOOST_UBLAS_INLINE
1184        const matrix_closure_type &data () const {
1185            return data_;
1186        }
1187        BOOST_UBLAS_INLINE
1188        matrix_closure_type &data () {
1189            return data_;
1190        }
1191
1192        // Element access
1193#ifndef BOOST_UBLAS_PROXY_CONST_MEMBER
1194        BOOST_UBLAS_INLINE
1195        const_reference operator () (size_type i, size_type j) const {
1196            BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1197            BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1198            // if (i == j)
1199            //     return type_traits<value_type>::real (data () (i, i));
1200            // else
1201            if (triangular_type::other (i, j))
1202                return data () (i, j);
1203            else
1204                return type_traits<value_type>::conj (data () (j, i));
1205        }
1206        BOOST_UBLAS_INLINE
1207        reference operator () (size_type i, size_type j) {
1208            BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1209            BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1210#ifndef BOOST_UBLAS_STRICT_HERMITIAN
1211            if (triangular_type::other (i, j))
1212                return data () (i, j);
1213            else {
1214                external_logic ().raise ();
1215                return conj_ = type_traits<value_type>::conj (data () (j, i));
1216            }
1217#else
1218            if (triangular_type::other (i, j))
1219                return reference (*this, i, j, data () (i, j));
1220            else
1221                return reference (*this, i, j, type_traits<value_type>::conj (data () (j, i)));
1222#endif
1223        }
1224        BOOST_UBLAS_INLINE
1225        true_reference insert_element (size_type i, size_type j, value_type t) {
1226            BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1227            BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1228            // if (i == j)
1229            //     data () (i, i) = type_traits<value_type>::real (t);
1230            // else
1231            if (triangular_type::other (i, j))
1232                return data () (i, j) = t;
1233            else
1234                return data () (j, i) = type_traits<value_type>::conj (t);
1235        }
1236#else
1237        BOOST_UBLAS_INLINE
1238        reference operator () (size_type i, size_type j) {
1239            BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1240            BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1241#ifndef BOOST_UBLAS_STRICT_HERMITIAN
1242            if (triangular_type::other (i, j))
1243                return data () (i, j);
1244            else {
1245                external_logic ().raise ();
1246                return conj_ = type_traits<value_type>::conj (data () (j, i));
1247            }
1248#else
1249            if (triangular_type::other (i, j))
1250                return reference (*this, i, j, data () (i, j));
1251            else
1252                return reference (*this, i, j, type_traits<value_type>::conj (data () (j, i)));
1253#endif
1254        }
1255        BOOST_UBLAS_INLINE
1256        true_reference insert_element (size_type i, size_type j, value_type t) {
1257            BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1258            BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1259            // if (i == j)
1260            //     data () (i, i) = type_traits<value_type>::real (t);
1261            // else
1262            if (triangular_type::other (i, j))
1263                return data () (i, j) = t;
1264            else
1265                return data () (j, i) = type_traits<value_type>::conj (t);
1266        }
1267#endif
1268
1269        // Assignment
1270        BOOST_UBLAS_INLINE
1271        hermitian_adaptor &operator = (const hermitian_adaptor &m) {
1272            matrix_assign<scalar_assign, triangular_type> (*this, m);
1273            return *this;
1274        }
1275        BOOST_UBLAS_INLINE
1276        hermitian_adaptor &assign_temporary (hermitian_adaptor &m) {
1277            *this = m;
1278            return *this;
1279        }
1280        template<class AE>
1281        BOOST_UBLAS_INLINE
1282        hermitian_adaptor &operator = (const matrix_expression<AE> &ae) {
1283            matrix_assign<scalar_assign, triangular_type> (*this, matrix<value_type> (ae));
1284            return *this;
1285        }
1286        template<class AE>
1287        BOOST_UBLAS_INLINE
1288        hermitian_adaptor &assign (const matrix_expression<AE> &ae) {
1289            matrix_assign<scalar_assign, triangular_type> (*this, ae);
1290            return *this;
1291        }
1292        template<class AE>
1293        BOOST_UBLAS_INLINE
1294        hermitian_adaptor& operator += (const matrix_expression<AE> &ae) {
1295            matrix_assign<scalar_assign, triangular_type> (*this, matrix<value_type> (*this + ae));
1296            return *this;
1297        }
1298        template<class AE>
1299        BOOST_UBLAS_INLINE
1300        hermitian_adaptor &plus_assign (const matrix_expression<AE> &ae) {
1301            matrix_assign<scalar_plus_assign, triangular_type> (*this, ae);
1302            return *this;
1303        }
1304        template<class AE>
1305        BOOST_UBLAS_INLINE
1306        hermitian_adaptor& operator -= (const matrix_expression<AE> &ae) {
1307            matrix_assign<scalar_assign, triangular_type> (*this, matrix<value_type> (*this - ae));
1308            return *this;
1309        }
1310        template<class AE>
1311        BOOST_UBLAS_INLINE
1312        hermitian_adaptor &minus_assign (const matrix_expression<AE> &ae) {
1313            matrix_assign<scalar_minus_assign, triangular_type> (*this, ae);
1314            return *this;
1315        }
1316        template<class AT>
1317        BOOST_UBLAS_INLINE
1318        hermitian_adaptor& operator *= (const AT &at) {
1319            // Multiplication is only allowed for real scalars,
1320            // otherwise the resulting matrix isn't hermitian.
1321            // Thanks to Peter Schmitteckert for spotting this.
1322            BOOST_UBLAS_CHECK (type_traits<value_type>::imag (at) == 0, non_real ());
1323            matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1324            return *this;
1325        }
1326        template<class AT>
1327        BOOST_UBLAS_INLINE
1328        hermitian_adaptor& operator /= (const AT &at) {
1329            // Multiplication is only allowed for real scalars,
1330            // otherwise the resulting matrix isn't hermitian.
1331            // Thanks to Peter Schmitteckert for spotting this.
1332            BOOST_UBLAS_CHECK (type_traits<value_type>::imag (at) == 0, non_real ());
1333            matrix_assign_scalar<scalar_divides_assign> (*this, at);
1334            return *this;
1335        }
1336
1337        // Closure comparison
1338        BOOST_UBLAS_INLINE
1339        bool same_closure (const hermitian_adaptor &ha) const {
1340            return (*this).data ().same_closure (ha.data ());
1341        }
1342
1343        // Swapping
1344        BOOST_UBLAS_INLINE
1345        void swap (hermitian_adaptor &m) {
1346            if (this != &m)
1347                matrix_swap<scalar_swap, triangular_type> (*this, m);
1348        }
1349        BOOST_UBLAS_INLINE
1350        friend void swap (hermitian_adaptor &m1, hermitian_adaptor &m2) {
1351            m1.swap (m2);
1352        }
1353
1354        // Iterator types
1355    private:
1356        // Use matrix iterator
1357        typedef typename M::const_iterator1 const_subiterator1_type;
1358        typedef typename boost::mpl::if_<boost::is_const<M>,
1359                                          typename M::const_iterator1,
1360                                          typename M::iterator1>::type subiterator1_type;
1361        typedef typename M::const_iterator2 const_subiterator2_type;
1362        typedef typename boost::mpl::if_<boost::is_const<M>,
1363                                          typename M::const_iterator2,
1364                                          typename M::iterator2>::type subiterator2_type;
1365
1366    public:
1367#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1368        typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
1369        typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
1370        typedef indexed_const_iterator1<self_type, dense_random_access_iterator_tag> const_iterator1;
1371        typedef indexed_const_iterator2<self_type, dense_random_access_iterator_tag> const_iterator2;
1372#else
1373        class const_iterator1;
1374        class iterator1;
1375        class const_iterator2;
1376        class iterator2;
1377#endif
1378        typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1379        typedef reverse_iterator_base1<iterator1> reverse_iterator1;
1380        typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1381        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
1382
1383        // Element lookup
1384        BOOST_UBLAS_INLINE
1385        const_iterator1 find1 (int rank, size_type i, size_type j) const {
1386            if (triangular_type::other (i, j)) {
1387                if (triangular_type::other (size1 (), j)) {
1388                    return const_iterator1 (*this, 0, 0,
1389                                            data ().find1 (rank, i, j), data ().find1 (rank, size1 (), j),
1390                                            data ().find2 (rank, size2 (), size1 ()), data ().find2 (rank, size2 (), size1 ()));
1391                } else {
1392                    return const_iterator1 (*this, 0, 1,
1393                                            data ().find1 (rank, i, j), data ().find1 (rank, j, j),
1394                                            data ().find2 (rank, j, j), data ().find2 (rank, j, size1 ()));
1395                }
1396            } else {
1397                if (triangular_type::other (size1 (), j)) {
1398                    return const_iterator1 (*this, 1, 0,
1399                                            data ().find1 (rank, j, j), data ().find1 (rank, size1 (), j),
1400                                            data ().find2 (rank, j, i), data ().find2 (rank, j, j));
1401                } else {
1402                    return const_iterator1 (*this, 1, 1,
1403                                            data ().find1 (rank, size1 (), size2 ()), data ().find1 (rank, size1 (), size2 ()),
1404                                            data ().find2 (rank, j, i), data ().find2 (rank, j, size1 ()));
1405                }
1406            }
1407        }
1408        BOOST_UBLAS_INLINE
1409        iterator1 find1 (int rank, size_type i, size_type j) {
1410            if (rank == 1)
1411                i = triangular_type::mutable_restrict1 (i, j);
1412            return iterator1 (*this, data ().find1 (rank, i, j));
1413        }
1414        BOOST_UBLAS_INLINE
1415        const_iterator2 find2 (int rank, size_type i, size_type j) const {
1416            if (triangular_type::other (i, j)) {
1417                if (triangular_type::other (i, size2 ())) {
1418                    return const_iterator2 (*this, 1, 1,
1419                                            data ().find1 (rank, size2 (), size1 ()), data ().find1 (rank, size2 (), size1 ()),
1420                                            data ().find2 (rank, i, j), data ().find2 (rank, i, size2 ()));
1421                } else {
1422                    return const_iterator2 (*this, 1, 0,
1423                                            data ().find1 (rank, i, i), data ().find1 (rank, size2 (), i),
1424                                            data ().find2 (rank, i, j), data ().find2 (rank, i, i));
1425                }
1426            } else {
1427                if (triangular_type::other (i, size2 ())) {
1428                    return const_iterator2 (*this, 0, 1,
1429                                            data ().find1 (rank, j, i), data ().find1 (rank, i, i),
1430                                            data ().find2 (rank, i, i), data ().find2 (rank, i, size2 ()));
1431                } else {
1432                    return const_iterator2 (*this, 0, 0,
1433                                            data ().find1 (rank, j, i), data ().find1 (rank, size2 (), i),
1434                                            data ().find2 (rank, size1 (), size2 ()), data ().find2 (rank, size2 (), size2 ()));
1435                }
1436            }
1437        }
1438        BOOST_UBLAS_INLINE
1439        iterator2 find2 (int rank, size_type i, size_type j) {
1440            if (rank == 1)
1441                j = triangular_type::mutable_restrict2 (i, j);
1442            return iterator2 (*this, data ().find2 (rank, i, j));
1443        }
1444
1445        // Iterators simply are indices.
1446
1447#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1448        class const_iterator1:
1449            public container_const_reference<hermitian_adaptor>,
1450            public random_access_iterator_base<typename iterator_restrict_traits<
1451                                                   typename const_subiterator1_type::iterator_category, dense_random_access_iterator_tag>::iterator_category,
1452                                               const_iterator1, value_type> {
1453        public:
1454            typedef typename const_subiterator1_type::value_type value_type;
1455            typedef typename const_subiterator1_type::difference_type difference_type;
1456            // FIXME no better way to not return the address of a temporary?
1457            // typedef typename const_subiterator1_type::reference reference;
1458            typedef typename const_subiterator1_type::value_type reference;
1459            typedef typename const_subiterator1_type::pointer pointer;
1460
1461            typedef const_iterator2 dual_iterator_type;
1462            typedef const_reverse_iterator2 dual_reverse_iterator_type;
1463
1464            // Construction and destruction
1465            BOOST_UBLAS_INLINE
1466            const_iterator1 ():
1467                container_const_reference<self_type> (),
1468                begin_ (-1), end_ (-1), current_ (-1),
1469                it1_begin_ (), it1_end_ (), it1_ (),
1470                it2_begin_ (), it2_end_ (), it2_ () {}
1471            BOOST_UBLAS_INLINE
1472            const_iterator1 (const self_type &m, int begin, int end,
1473                             const const_subiterator1_type &it1_begin, const const_subiterator1_type &it1_end,
1474                             const const_subiterator2_type &it2_begin, const const_subiterator2_type &it2_end):
1475                container_const_reference<self_type> (m),
1476                begin_ (begin), end_ (end), current_ (begin),
1477                it1_begin_ (it1_begin), it1_end_ (it1_end), it1_ (it1_begin_),
1478                it2_begin_ (it2_begin), it2_end_ (it2_end), it2_ (it2_begin_) {
1479                if (current_ == 0 && it1_ == it1_end_)
1480                    current_ = 1;
1481                if (current_ == 1 && it2_ == it2_end_)
1482                    current_ = 0;
1483                if ((current_ == 0 && it1_ == it1_end_) ||
1484                    (current_ == 1 && it2_ == it2_end_))
1485                    current_ = end_;
1486                BOOST_UBLAS_CHECK (current_ == end_ ||
1487                                   (current_ == 0 && it1_ != it1_end_) ||
1488                                   (current_ == 1 && it2_ != it2_end_), internal_logic ());
1489            }
1490            // FIXME cannot compile
1491            //  iterator1 does not have these members!
1492            BOOST_UBLAS_INLINE
1493            const_iterator1 (const iterator1 &it):
1494                container_const_reference<self_type> (it ()),
1495                begin_ (it.begin_), end_ (it.end_), current_ (it.current_),
1496                it1_begin_ (it.it1_begin_), it1_end_ (it.it1_end_), it1_ (it.it1_),
1497                it2_begin_ (it.it2_begin_), it2_end_ (it.it2_end_), it2_ (it.it2_) {
1498                BOOST_UBLAS_CHECK (current_ == end_ ||
1499                                   (current_ == 0 && it1_ != it1_end_) ||
1500                                   (current_ == 1 && it2_ != it2_end_), internal_logic ());
1501            }
1502
1503            // Arithmetic
1504            BOOST_UBLAS_INLINE
1505            const_iterator1 &operator ++ () {
1506                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1507                if (current_ == 0) {
1508                    BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1509                    ++ it1_;
1510                    if (it1_ == it1_end_ && end_ == 1) {
1511                        it2_ = it2_begin_;
1512                        current_ = 1;
1513                    }
1514                } else /* if (current_ == 1) */ {
1515                    BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1516                    ++ it2_;
1517                    if (it2_ == it2_end_ && end_ == 0) {
1518                        it1_ = it1_begin_;
1519                        current_ = 0;
1520                    }
1521                }
1522                return *this;
1523            }
1524            BOOST_UBLAS_INLINE
1525            const_iterator1 &operator -- () {
1526                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1527                if (current_ == 0) {
1528                    if (it1_ == it1_begin_ && begin_ == 1) {
1529                        it2_ = it2_end_;
1530                        BOOST_UBLAS_CHECK (it2_ != it2_begin_, internal_logic ());
1531                        -- it2_;
1532                        current_ = 1;
1533                    } else {
1534                        -- it1_;
1535                    }
1536                } else /* if (current_ == 1) */ {
1537                    if (it2_ == it2_begin_ && begin_ == 0) {
1538                        it1_ = it1_end_;
1539                        BOOST_UBLAS_CHECK (it1_ != it1_begin_, internal_logic ());
1540                        -- it1_;
1541                        current_ = 0;
1542                    } else {
1543                        -- it2_;
1544                    }
1545                }
1546                return *this;
1547            }
1548            BOOST_UBLAS_INLINE
1549            const_iterator1 &operator += (difference_type n) {
1550                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1551                if (current_ == 0) {
1552                    size_type d = (std::min) (n, it1_end_ - it1_);
1553                    it1_ += d;
1554                    n -= d;
1555                    if (n > 0 || (end_ == 1 && it1_ == it1_end_)) {
1556                        BOOST_UBLAS_CHECK (end_ == 1, external_logic ());
1557                        d = (std::min) (n, it2_end_ - it2_begin_);
1558                        it2_ = it2_begin_ + d;
1559                        n -= d;
1560                        current_ = 1;
1561                    }
1562                } else /* if (current_ == 1) */ {
1563                    size_type d = (std::min) (n, it2_end_ - it2_);
1564                    it2_ += d;
1565                    n -= d;
1566                    if (n > 0 || (end_ == 0 && it2_ == it2_end_)) {
1567                        BOOST_UBLAS_CHECK (end_ == 0, external_logic ());
1568                        d = (std::min) (n, it1_end_ - it1_begin_);
1569                        it1_ = it1_begin_ + d;
1570                        n -= d;
1571                        current_ = 0;
1572                    }
1573                }
1574                BOOST_UBLAS_CHECK (n == 0, external_logic ());
1575                return *this;
1576            }
1577            BOOST_UBLAS_INLINE
1578            const_iterator1 &operator -= (difference_type n) {
1579                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1580                if (current_ == 0) {
1581                    size_type d = (std::min) (n, it1_ - it1_begin_);
1582                    it1_ -= d;
1583                    n -= d;
1584                    if (n > 0) {
1585                        BOOST_UBLAS_CHECK (end_ == 1, external_logic ());
1586                        d = (std::min) (n, it2_end_ - it2_begin_);
1587                        it2_ = it2_end_ - d;
1588                        n -= d;
1589                        current_ = 1;
1590                    }
1591                } else /* if (current_ == 1) */ {
1592                    size_type d = (std::min) (n, it2_ - it2_begin_);
1593                    it2_ -= d;
1594                    n -= d;
1595                    if (n > 0) {
1596                        BOOST_UBLAS_CHECK (end_ == 0, external_logic ());
1597                        d = (std::min) (n, it1_end_ - it1_begin_);
1598                        it1_ = it1_end_ - d;
1599                        n -= d;
1600                        current_ = 0;
1601                    }
1602                }
1603                BOOST_UBLAS_CHECK (n == 0, external_logic ());
1604                return *this;
1605            }
1606            BOOST_UBLAS_INLINE
1607            difference_type operator - (const const_iterator1 &it) const {
1608                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1609                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1610                BOOST_UBLAS_CHECK (it.current_ == 0 || it.current_ == 1, internal_logic ());
1611                BOOST_UBLAS_CHECK (/* begin_ == it.begin_ && */ end_ == it.end_, internal_logic ());
1612                if (current_ == 0 && it.current_ == 0) {
1613                    return it1_ - it.it1_;
1614                } else if (current_ == 0 && it.current_ == 1) {
1615                    if (end_ == 1 && it.end_ == 1) {
1616                        return (it1_ - it.it1_end_) + (it.it2_begin_ - it.it2_);
1617                    } else /* if (end_ == 0 && it.end_ == 0) */ {
1618                        return (it1_ - it.it1_begin_) + (it.it2_end_ - it.it2_);
1619                    }
1620
1621                } else if (current_ == 1 && it.current_ == 0) {
1622                    if (end_ == 1 && it.end_ == 1) {
1623                        return (it2_ - it.it2_begin_) + (it.it1_end_ - it.it1_);
1624                    } else /* if (end_ == 0 && it.end_ == 0) */ {
1625                        return (it2_ - it.it2_end_) + (it.it1_begin_ - it.it1_);
1626                    }
1627                } else /* if (current_ == 1 && it.current_ == 1) */ {
1628                    return it2_ - it.it2_;
1629                }
1630            }
1631
1632            // Dereference
1633            BOOST_UBLAS_INLINE
1634            const_reference operator * () const {
1635                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1636                if (current_ == 0) {
1637                    BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1638                    if (triangular_type::other (index1 (), index2 ()))
1639                        return *it1_;
1640                    else
1641                        return type_traits<value_type>::conj (*it1_);
1642                } else /* if (current_ == 1) */ {
1643                    BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1644                    if (triangular_type::other (index1 (), index2 ()))
1645                        return *it2_;
1646                    else
1647                        return type_traits<value_type>::conj (*it2_);
1648                }
1649            }
1650
1651#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1652            BOOST_UBLAS_INLINE
1653#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1654            typename self_type::
1655#endif
1656            const_iterator2 begin () const {
1657                return (*this) ().find2 (1, index1 (), 0);
1658            }
1659            BOOST_UBLAS_INLINE
1660#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1661            typename self_type::
1662#endif
1663            const_iterator2 end () const {
1664                return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1665            }
1666            BOOST_UBLAS_INLINE
1667#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1668            typename self_type::
1669#endif
1670            const_reverse_iterator2 rbegin () const {
1671                return const_reverse_iterator2 (end ());
1672            }
1673            BOOST_UBLAS_INLINE
1674#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1675            typename self_type::
1676#endif
1677            const_reverse_iterator2 rend () const {
1678                return const_reverse_iterator2 (begin ());
1679            }
1680#endif
1681
1682            // Indices
1683            BOOST_UBLAS_INLINE
1684            size_type index1 () const {
1685                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1686                if (current_ == 0) {
1687                    BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1688                    return it1_.index1 ();
1689                } else /* if (current_ == 1) */ {
1690                    BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1691                    return it2_.index2 ();
1692                }
1693            }
1694            BOOST_UBLAS_INLINE
1695            size_type index2 () const {
1696                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1697                if (current_ == 0) {
1698                    BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1699                    return it1_.index2 ();
1700                } else /* if (current_ == 1) */ {
1701                    BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1702                    return it2_.index1 ();
1703                }
1704            }
1705
1706            // Assignment
1707            BOOST_UBLAS_INLINE
1708            const_iterator1 &operator = (const const_iterator1 &it) {
1709                container_const_reference<self_type>::assign (&it ());
1710                begin_ = it.begin_;
1711                end_ = it.end_;
1712                current_ = it.current_;
1713                it1_begin_ = it.it1_begin_;
1714                it1_end_ = it.it1_end_;
1715                it1_ = it.it1_;
1716                it2_begin_ = it.it2_begin_;
1717                it2_end_ = it.it2_end_;
1718                it2_ = it.it2_;
1719                return *this;
1720            }
1721
1722            // Comparison
1723            BOOST_UBLAS_INLINE
1724            bool operator == (const const_iterator1 &it) const {
1725                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1726                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1727                BOOST_UBLAS_CHECK (it.current_ == 0 || it.current_ == 1, internal_logic ());
1728                BOOST_UBLAS_CHECK (/* begin_ == it.begin_ && */ end_ == it.end_, internal_logic ());
1729                return (current_ == 0 && it.current_ == 0 && it1_ == it.it1_) ||
1730                       (current_ == 1 && it.current_ == 1 && it2_ == it.it2_);
1731            }
1732            BOOST_UBLAS_INLINE
1733            bool operator < (const const_iterator1 &it) const {
1734                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1735                return it - *this > 0;
1736            }
1737
1738        private:
1739            int begin_;
1740            int end_;
1741            int current_;
1742            const_subiterator1_type it1_begin_;
1743            const_subiterator1_type it1_end_;
1744            const_subiterator1_type it1_;
1745            const_subiterator2_type it2_begin_;
1746            const_subiterator2_type it2_end_;
1747            const_subiterator2_type it2_;
1748        };
1749#endif
1750
1751        BOOST_UBLAS_INLINE
1752        const_iterator1 begin1 () const {
1753            return find1 (0, 0, 0);
1754        }
1755        BOOST_UBLAS_INLINE
1756        const_iterator1 end1 () const {
1757            return find1 (0, size1 (), 0);
1758        }
1759
1760#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1761        class iterator1:
1762            public container_reference<hermitian_adaptor>,
1763            public random_access_iterator_base<typename iterator_restrict_traits<
1764                                                   typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1765                                               iterator1, value_type> {
1766        public:
1767            typedef typename subiterator1_type::value_type value_type;
1768            typedef typename subiterator1_type::difference_type difference_type;
1769            typedef typename subiterator1_type::reference reference;
1770            typedef typename subiterator1_type::pointer pointer;
1771
1772            typedef iterator2 dual_iterator_type;
1773            typedef reverse_iterator2 dual_reverse_iterator_type;
1774
1775            // Construction and destruction
1776            BOOST_UBLAS_INLINE
1777            iterator1 ():
1778                container_reference<self_type> (), it1_ () {}
1779            BOOST_UBLAS_INLINE
1780            iterator1 (self_type &m, const subiterator1_type &it1):
1781                container_reference<self_type> (m), it1_ (it1) {}
1782
1783            // Arithmetic
1784            BOOST_UBLAS_INLINE
1785            iterator1 &operator ++ () {
1786                ++ it1_;
1787                return *this;
1788            }
1789            BOOST_UBLAS_INLINE
1790            iterator1 &operator -- () {
1791                -- it1_;
1792                return *this;
1793            }
1794            BOOST_UBLAS_INLINE
1795            iterator1 &operator += (difference_type n) {
1796                it1_ += n;
1797                return *this;
1798            }
1799            BOOST_UBLAS_INLINE
1800            iterator1 &operator -= (difference_type n) {
1801                it1_ -= n;
1802                return *this;
1803            }
1804            BOOST_UBLAS_INLINE
1805            difference_type operator - (const iterator1 &it) const {
1806                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1807                return it1_ - it.it1_;
1808            }
1809
1810            // Dereference
1811            BOOST_UBLAS_INLINE
1812            reference operator * () const {
1813                return *it1_;
1814            }
1815
1816#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1817            BOOST_UBLAS_INLINE
1818#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1819            typename self_type::
1820#endif
1821            iterator2 begin () const {
1822                return (*this) ().find2 (1, index1 (), 0);
1823            }
1824            BOOST_UBLAS_INLINE
1825#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1826            typename self_type::
1827#endif
1828            iterator2 end () const {
1829                return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1830            }
1831            BOOST_UBLAS_INLINE
1832#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1833            typename self_type::
1834#endif
1835            reverse_iterator2 rbegin () const {
1836                return reverse_iterator2 (end ());
1837            }
1838            BOOST_UBLAS_INLINE
1839#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1840            typename self_type::
1841#endif
1842            reverse_iterator2 rend () const {
1843                return reverse_iterator2 (begin ());
1844            }
1845#endif
1846
1847            // Indices
1848            BOOST_UBLAS_INLINE
1849            size_type index1 () const {
1850                return it1_.index1 ();
1851            }
1852            BOOST_UBLAS_INLINE
1853            size_type index2 () const {
1854                return it1_.index2 ();
1855            }
1856
1857            // Assignment
1858            BOOST_UBLAS_INLINE
1859            iterator1 &operator = (const iterator1 &it) {
1860                container_reference<self_type>::assign (&it ());
1861                it1_ = it.it1_;
1862                return *this;
1863            }
1864
1865            // Comparison
1866            BOOST_UBLAS_INLINE
1867            bool operator == (const iterator1 &it) const {
1868                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1869                return it1_ == it.it1_;
1870            }
1871            BOOST_UBLAS_INLINE
1872            bool operator < (const iterator1 &it) const {
1873                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1874                return it1_ < it.it1_;
1875            }
1876
1877        private:
1878            subiterator1_type it1_;
1879
1880            friend class const_iterator1;
1881        };
1882#endif
1883
1884        BOOST_UBLAS_INLINE
1885        iterator1 begin1 () {
1886            return find1 (0, 0, 0);
1887        }
1888        BOOST_UBLAS_INLINE
1889        iterator1 end1 () {
1890            return find1 (0, size1 (), 0);
1891        }
1892
1893#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1894        class const_iterator2:
1895            public container_const_reference<hermitian_adaptor>,
1896            public random_access_iterator_base<typename iterator_restrict_traits<
1897            typename const_subiterator2_type::iterator_category, dense_random_access_iterator_tag>::iterator_category,
1898                                               const_iterator2, value_type> {
1899        public:
1900            typedef typename const_subiterator2_type::value_type value_type;
1901            typedef typename const_subiterator2_type::difference_type difference_type;
1902            // FIXME no better way to not return the address of a temporary?
1903            // typedef typename const_subiterator2_type::reference reference;
1904            typedef typename const_subiterator2_type::value_type reference;
1905            typedef typename const_subiterator2_type::pointer pointer;
1906
1907            typedef const_iterator1 dual_iterator_type;
1908            typedef const_reverse_iterator1 dual_reverse_iterator_type;
1909
1910            // Construction and destruction
1911            BOOST_UBLAS_INLINE
1912            const_iterator2 ():
1913                container_const_reference<self_type> (),
1914                begin_ (-1), end_ (-1), current_ (-1),
1915                it1_begin_ (), it1_end_ (), it1_ (),
1916                it2_begin_ (), it2_end_ (), it2_ () {}
1917            BOOST_UBLAS_INLINE
1918            const_iterator2 (const self_type &m, int begin, int end,
1919                             const const_subiterator1_type &it1_begin, const const_subiterator1_type &it1_end,
1920                             const const_subiterator2_type &it2_begin, const const_subiterator2_type &it2_end):
1921                container_const_reference<self_type> (m),
1922                begin_ (begin), end_ (end), current_ (begin),
1923                it1_begin_ (it1_begin), it1_end_ (it1_end), it1_ (it1_begin_),
1924                it2_begin_ (it2_begin), it2_end_ (it2_end), it2_ (it2_begin_) {
1925                if (current_ == 0 && it1_ == it1_end_)
1926                    current_ = 1;
1927                if (current_ == 1 && it2_ == it2_end_)
1928                    current_ = 0;
1929                if ((current_ == 0 && it1_ == it1_end_) ||
1930                    (current_ == 1 && it2_ == it2_end_))
1931                    current_ = end_;
1932                BOOST_UBLAS_CHECK (current_ == end_ ||
1933                                   (current_ == 0 && it1_ != it1_end_) ||
1934                                   (current_ == 1 && it2_ != it2_end_), internal_logic ());
1935            }
1936            // FIXME cannot compiler
1937            //  iterator2 does not have these members!
1938            BOOST_UBLAS_INLINE
1939            const_iterator2 (const iterator2 &it):
1940                container_const_reference<self_type> (it ()),
1941                begin_ (it.begin_), end_ (it.end_), current_ (it.current_),
1942                it1_begin_ (it.it1_begin_), it1_end_ (it.it1_end_), it1_ (it.it1_),
1943                it2_begin_ (it.it2_begin_), it2_end_ (it.it2_end_), it2_ (it.it2_) {
1944                BOOST_UBLAS_CHECK (current_ == end_ ||
1945                                   (current_ == 0 && it1_ != it1_end_) ||
1946                                   (current_ == 1 && it2_ != it2_end_), internal_logic ());
1947            }
1948
1949            // Arithmetic
1950            BOOST_UBLAS_INLINE
1951            const_iterator2 &operator ++ () {
1952                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1953                if (current_ == 0) {
1954                    BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1955                    ++ it1_;
1956                    if (it1_ == it1_end_ && end_ == 1) {
1957                        it2_ = it2_begin_;
1958                        current_ = 1;
1959                    }
1960                } else /* if (current_ == 1) */ {
1961                    BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1962                    ++ it2_;
1963                    if (it2_ == it2_end_ && end_ == 0) {
1964                        it1_ = it1_begin_;
1965                        current_ = 0;
1966                    }
1967                }
1968                return *this;
1969            }
1970            BOOST_UBLAS_INLINE
1971            const_iterator2 &operator -- () {
1972                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1973                if (current_ == 0) {
1974                    if (it1_ == it1_begin_ && begin_ == 1) {
1975                        it2_ = it2_end_;
1976                        BOOST_UBLAS_CHECK (it2_ != it2_begin_, internal_logic ());
1977                        -- it2_;
1978                        current_ = 1;
1979                    } else {
1980                        -- it1_;
1981                    }
1982                } else /* if (current_ == 1) */ {
1983                    if (it2_ == it2_begin_ && begin_ == 0) {
1984                        it1_ = it1_end_;
1985                        BOOST_UBLAS_CHECK (it1_ != it1_begin_, internal_logic ());
1986                        -- it1_;
1987                        current_ = 0;
1988                    } else {
1989                        -- it2_;
1990                    }
1991                }
1992                return *this;
1993            }
1994            BOOST_UBLAS_INLINE
1995            const_iterator2 &operator += (difference_type n) {
1996                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1997                if (current_ == 0) {
1998                    size_type d = (std::min) (n, it1_end_ - it1_);
1999                    it1_ += d;
2000                    n -= d;
2001                    if (n > 0 || (end_ == 1 && it1_ == it1_end_)) {
2002                        BOOST_UBLAS_CHECK (end_ == 1, external_logic ());
2003                        d = (std::min) (n, it2_end_ - it2_begin_);
2004                        it2_ = it2_begin_ + d;
2005                        n -= d;
2006                        current_ = 1;
2007                    }
2008                } else /* if (current_ == 1) */ {
2009                    size_type d = (std::min) (n, it2_end_ - it2_);
2010                    it2_ += d;
2011                    n -= d;
2012                    if (n > 0 || (end_ == 0 && it2_ == it2_end_)) {
2013                        BOOST_UBLAS_CHECK (end_ == 0, external_logic ());
2014                        d = (std::min) (n, it1_end_ - it1_begin_);
2015                        it1_ = it1_begin_ + d;
2016                        n -= d;
2017                        current_ = 0;
2018                    }
2019                }
2020                BOOST_UBLAS_CHECK (n == 0, external_logic ());
2021                return *this;
2022            }
2023            BOOST_UBLAS_INLINE
2024            const_iterator2 &operator -= (difference_type n) {
2025                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2026                if (current_ == 0) {
2027                    size_type d = (std::min) (n, it1_ - it1_begin_);
2028                    it1_ -= d;
2029                    n -= d;
2030                    if (n > 0) {
2031                        BOOST_UBLAS_CHECK (end_ == 1, external_logic ());
2032                        d = (std::min) (n, it2_end_ - it2_begin_);
2033                        it2_ = it2_end_ - d;
2034                        n -= d;
2035                        current_ = 1;
2036                    }
2037                } else /* if (current_ == 1) */ {
2038                    size_type d = (std::min) (n, it2_ - it2_begin_);
2039                    it2_ -= d;
2040                    n -= d;
2041                    if (n > 0) {
2042                        BOOST_UBLAS_CHECK (end_ == 0, external_logic ());
2043                        d = (std::min) (n, it1_end_ - it1_begin_);
2044                        it1_ = it1_end_ - d;
2045                        n -= d;
2046                        current_ = 0;
2047                    }
2048                }
2049                BOOST_UBLAS_CHECK (n == 0, external_logic ());
2050                return *this;
2051            }
2052            BOOST_UBLAS_INLINE
2053            difference_type operator - (const const_iterator2 &it) const {
2054                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2055                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2056                BOOST_UBLAS_CHECK (it.current_ == 0 || it.current_ == 1, internal_logic ());
2057                BOOST_UBLAS_CHECK (/* begin_ == it.begin_ && */ end_ == it.end_, internal_logic ());
2058                if (current_ == 0 && it.current_ == 0) {
2059                    return it1_ - it.it1_;
2060                } else if (current_ == 0 && it.current_ == 1) {
2061                    if (end_ == 1 && it.end_ == 1) {
2062                        return (it1_ - it.it1_end_) + (it.it2_begin_ - it.it2_);
2063                    } else /* if (end_ == 0 && it.end_ == 0) */ {
2064                        return (it1_ - it.it1_begin_) + (it.it2_end_ - it.it2_);
2065                    }
2066
2067                } else if (current_ == 1 && it.current_ == 0) {
2068                    if (end_ == 1 && it.end_ == 1) {
2069                        return (it2_ - it.it2_begin_) + (it.it1_end_ - it.it1_);
2070                    } else /* if (end_ == 0 && it.end_ == 0) */ {
2071                        return (it2_ - it.it2_end_) + (it.it1_begin_ - it.it1_);
2072                    }
2073                } else /* if (current_ == 1 && it.current_ == 1) */ {
2074                    return it2_ - it.it2_;
2075                }
2076            }
2077
2078            // Dereference
2079            BOOST_UBLAS_INLINE
2080            const_reference operator * () const {
2081                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2082                if (current_ == 0) {
2083                    BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
2084                    if (triangular_type::other (index1 (), index2 ()))
2085                        return *it1_;
2086                    else
2087                        return type_traits<value_type>::conj (*it1_);
2088                } else /* if (current_ == 1) */ {
2089                    BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
2090                    if (triangular_type::other (index1 (), index2 ()))
2091                        return *it2_;
2092                    else
2093                        return type_traits<value_type>::conj (*it2_);
2094                }
2095            }
2096
2097#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2098            BOOST_UBLAS_INLINE
2099#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2100            typename self_type::
2101#endif
2102            const_iterator1 begin () const {
2103                return (*this) ().find1 (1, 0, index2 ());
2104            }
2105            BOOST_UBLAS_INLINE
2106#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2107            typename self_type::
2108#endif
2109            const_iterator1 end () const {
2110                return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2111            }
2112            BOOST_UBLAS_INLINE
2113#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2114            typename self_type::
2115#endif
2116            const_reverse_iterator1 rbegin () const {
2117                return const_reverse_iterator1 (end ());
2118            }
2119            BOOST_UBLAS_INLINE
2120#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2121            typename self_type::
2122#endif
2123            const_reverse_iterator1 rend () const {
2124                return const_reverse_iterator1 (begin ());
2125            }
2126#endif
2127
2128            // Indices
2129            BOOST_UBLAS_INLINE
2130            size_type index1 () const {
2131                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2132                if (current_ == 0) {
2133                    BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
2134                    return it1_.index2 ();
2135                } else /* if (current_ == 1) */ {
2136                    BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
2137                    return it2_.index1 ();
2138                }
2139            }
2140            BOOST_UBLAS_INLINE
2141            size_type index2 () const {
2142                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2143                if (current_ == 0) {
2144                    BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
2145                    return it1_.index1 ();
2146                } else /* if (current_ == 1) */ {
2147                    BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
2148                    return it2_.index2 ();
2149                }
2150            }
2151
2152            // Assignment
2153            BOOST_UBLAS_INLINE
2154            const_iterator2 &operator = (const const_iterator2 &it) {
2155                container_const_reference<self_type>::assign (&it ());
2156                begin_ = it.begin_;
2157                end_ = it.end_;
2158                current_ = it.current_;
2159                it1_begin_ = it.it1_begin_;
2160                it1_end_ = it.it1_end_;
2161                it1_ = it.it1_;
2162                it2_begin_ = it.it2_begin_;
2163                it2_end_ = it.it2_end_;
2164                it2_ = it.it2_;
2165                return *this;
2166            }
2167
2168            // Comparison
2169            BOOST_UBLAS_INLINE
2170            bool operator == (const const_iterator2 &it) const {
2171                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2172                BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2173                BOOST_UBLAS_CHECK (it.current_ == 0 || it.current_ == 1, internal_logic ());
2174                BOOST_UBLAS_CHECK (/* begin_ == it.begin_ && */ end_ == it.end_, internal_logic ());
2175                return (current_ == 0 && it.current_ == 0 && it1_ == it.it1_) ||
2176                       (current_ == 1 && it.current_ == 1 && it2_ == it.it2_);
2177            }
2178            BOOST_UBLAS_INLINE
2179            bool operator < (const const_iterator2 &it) const {
2180                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2181                return it - *this > 0;
2182            }
2183
2184        private:
2185            int begin_;
2186            int end_;
2187            int current_;
2188            const_subiterator1_type it1_begin_;
2189            const_subiterator1_type it1_end_;
2190            const_subiterator1_type it1_;
2191            const_subiterator2_type it2_begin_;
2192            const_subiterator2_type it2_end_;
2193            const_subiterator2_type it2_;
2194        };
2195#endif
2196
2197        BOOST_UBLAS_INLINE
2198        const_iterator2 begin2 () const {
2199            return find2 (0, 0, 0);
2200        }
2201        BOOST_UBLAS_INLINE
2202        const_iterator2 end2 () const {
2203            return find2 (0, 0, size2 ());
2204        }
2205
2206#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2207        class iterator2:
2208            public container_reference<hermitian_adaptor>,
2209            public random_access_iterator_base<typename iterator_restrict_traits<
2210                                                   typename subiterator2_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
2211                                               iterator2, value_type> {
2212        public:
2213            typedef typename subiterator2_type::value_type value_type;
2214            typedef typename subiterator2_type::difference_type difference_type;
2215            typedef typename subiterator2_type::reference reference;
2216            typedef typename subiterator2_type::pointer pointer;
2217
2218            typedef iterator1 dual_iterator_type;
2219            typedef reverse_iterator1 dual_reverse_iterator_type;
2220
2221            // Construction and destruction
2222            BOOST_UBLAS_INLINE
2223            iterator2 ():
2224                container_reference<self_type> (), it2_ () {}
2225            BOOST_UBLAS_INLINE
2226            iterator2 (self_type &m, const subiterator2_type &it2):
2227                container_reference<self_type> (m), it2_ (it2) {}
2228
2229            // Arithmetic
2230            BOOST_UBLAS_INLINE
2231            iterator2 &operator ++ () {
2232                ++ it2_;
2233                return *this;
2234            }
2235            BOOST_UBLAS_INLINE
2236            iterator2 &operator -- () {
2237                -- it2_;
2238                return *this;
2239            }
2240            BOOST_UBLAS_INLINE
2241            iterator2 &operator += (difference_type n) {
2242                it2_ += n;
2243                return *this;
2244            }
2245            BOOST_UBLAS_INLINE
2246            iterator2 &operator -= (difference_type n) {
2247                it2_ -= n;
2248                return *this;
2249            }
2250            BOOST_UBLAS_INLINE
2251            difference_type operator - (const iterator2 &it) const {
2252                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2253                return it2_ - it.it2_;
2254            }
2255
2256            // Dereference
2257            BOOST_UBLAS_INLINE
2258            reference operator * () const {
2259                return *it2_;
2260            }
2261
2262#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2263            BOOST_UBLAS_INLINE
2264#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2265            typename self_type::
2266#endif
2267            iterator1 begin () const {
2268                return (*this) ().find1 (1, 0, index2 ());
2269            }
2270            BOOST_UBLAS_INLINE
2271#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2272            typename self_type::
2273#endif
2274            iterator1 end () const {
2275                return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2276            }
2277            BOOST_UBLAS_INLINE
2278#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2279            typename self_type::
2280#endif
2281            reverse_iterator1 rbegin () const {
2282                return reverse_iterator1 (end ());
2283            }
2284            BOOST_UBLAS_INLINE
2285#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2286            typename self_type::
2287#endif
2288            reverse_iterator1 rend () const {
2289                return reverse_iterator1 (begin ());
2290            }
2291#endif
2292
2293            // Indices
2294            BOOST_UBLAS_INLINE
2295            size_type index1 () const {
2296                return it2_.index1 ();
2297            }
2298            BOOST_UBLAS_INLINE
2299            size_type index2 () const {
2300                return it2_.index2 ();
2301            }
2302
2303            // Assignment
2304            BOOST_UBLAS_INLINE
2305            iterator2 &operator = (const iterator2 &it) {
2306                container_reference<self_type>::assign (&it ());
2307                it2_ = it.it2_;
2308                return *this;
2309            }
2310
2311            // Comparison
2312            BOOST_UBLAS_INLINE
2313            bool operator == (const iterator2 &it) const {
2314                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2315                return it2_ == it.it2_;
2316            }
2317            BOOST_UBLAS_INLINE
2318            bool operator < (const iterator2 &it) const {
2319                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2320                return it2_ < it.it2_;
2321            }
2322
2323        private:
2324            subiterator2_type it2_;
2325
2326            friend class const_iterator2;
2327        };
2328#endif
2329
2330        BOOST_UBLAS_INLINE
2331        iterator2 begin2 () {
2332            return find2 (0, 0, 0);
2333        }
2334        BOOST_UBLAS_INLINE
2335        iterator2 end2 () {
2336            return find2 (0, 0, size2 ());
2337        }
2338
2339        // Reverse iterators
2340
2341        BOOST_UBLAS_INLINE
2342        const_reverse_iterator1 rbegin1 () const {
2343            return const_reverse_iterator1 (end1 ());
2344        }
2345        BOOST_UBLAS_INLINE
2346        const_reverse_iterator1 rend1 () const {
2347            return const_reverse_iterator1 (begin1 ());
2348        }
2349
2350        BOOST_UBLAS_INLINE
2351        reverse_iterator1 rbegin1 () {
2352            return reverse_iterator1 (end1 ());
2353        }
2354        BOOST_UBLAS_INLINE
2355        reverse_iterator1 rend1 () {
2356            return reverse_iterator1 (begin1 ());
2357        }
2358
2359        BOOST_UBLAS_INLINE
2360        const_reverse_iterator2 rbegin2 () const {
2361            return const_reverse_iterator2 (end2 ());
2362        }
2363        BOOST_UBLAS_INLINE
2364        const_reverse_iterator2 rend2 () const {
2365            return const_reverse_iterator2 (begin2 ());
2366        }
2367
2368        BOOST_UBLAS_INLINE
2369        reverse_iterator2 rbegin2 () {
2370            return reverse_iterator2 (end2 ());
2371        }
2372        BOOST_UBLAS_INLINE
2373        reverse_iterator2 rend2 () {
2374            return reverse_iterator2 (begin2 ());
2375        }
2376
2377    private:
2378        matrix_closure_type data_;
2379        static value_type conj_;
2380    };
2381
2382    template<class M, class TRI>
2383    typename hermitian_adaptor<M, TRI>::value_type hermitian_adaptor<M, TRI>::conj_;
2384
2385    // Specialization for temporary_traits
2386    template <class M, class TRI>
2387    struct vector_temporary_traits< hermitian_adaptor<M, TRI> >
2388    : vector_temporary_traits< M > {} ;
2389
2390    template <class M, class TRI>
2391    struct matrix_temporary_traits< hermitian_adaptor<M, TRI> >
2392    : matrix_temporary_traits< M > {} ;
2393
2394}}}
2395
2396#endif
Note: See TracBrowser for help on using the repository browser.