Update contrib.
2 // Copyright (c) 2000-2002
3 // Joerg Walter, Mathias Koch
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.
13 // The authors gratefully acknowledge the support of
14 // GeNeSys mbH & Co. KG in producing this work.
17 #ifndef BOOST_UBLAS_HERMITIAN_H
18 #define BOOST_UBLAS_HERMITIAN_H
20 #include <boost/numeric/ublas/matrix.hpp>
21 #include <boost/numeric/ublas/detail/temporary.hpp>
23 // Iterators based on ideas of Jeremy Siek
24 // Hermitian matrices are square. Thanks to Peter Schmitteckert for spotting this.
26 namespace boost { namespace numeric { namespace ublas {
29 bool is_hermitian (const M &m) {
30 typedef typename M::size_type size_type;
32 if (m.size1 () != m.size2 ())
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)))
44 #ifdef BOOST_UBLAS_STRICT_HERMITIAN
47 class hermitian_matrix_element:
48 public container_reference<M> {
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;
57 // Construction and destruction
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) {}
62 hermitian_matrix_element (const hermitian_matrix_element &p):
63 container_reference<matrix_type> (p), i_ (p.i_), d_ (p.d_), dirty_ (p.dirty_) {}
65 ~hermitian_matrix_element () {
67 (*this) ().insert_element (i_, j_, d_);
72 hermitian_matrix_element &operator = (const hermitian_matrix_element &p) {
73 // Overide the implict copy assignment
80 hermitian_matrix_element &operator = (const D &d) {
87 hermitian_matrix_element &operator += (const D &d) {
94 hermitian_matrix_element &operator -= (const D &d) {
101 hermitian_matrix_element &operator *= (const D &d) {
108 hermitian_matrix_element &operator /= (const D &d) {
117 bool operator == (const D &d) const {
122 bool operator != (const D &d) const {
128 operator const_reference () const {
134 void swap (hermitian_matrix_element p) {
138 std::swap (d_, p.d_);
142 friend void swap (hermitian_matrix_element p1, hermitian_matrix_element p2) {
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;
163 static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
164 static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
168 real_type real (const_reference t) {
169 return type_traits<element_type>::real (t);
173 real_type imag (const_reference t) {
174 return type_traits<element_type>::imag (t);
178 value_type conj (const_reference t) {
179 return type_traits<element_type>::conj (t);
184 real_type type_abs (const_reference t) {
185 return type_traits<element_type>::type_abs (t);
189 value_type type_sqrt (const_reference t) {
190 return type_traits<element_type>::type_sqrt (t);
195 real_type norm_1 (const_reference t) {
196 return type_traits<element_type>::norm_1 (t);
200 real_type norm_2 (const_reference t) {
201 return type_traits<element_type>::norm_2 (t);
205 real_type norm_inf (const_reference t) {
206 return type_traits<element_type>::norm_inf (t);
211 bool equals (const_reference t1, const_reference t2) {
212 return type_traits<element_type>::equals (t1, t2);
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;
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;
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;
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> > {
237 typedef T &true_reference;
239 typedef TRI triangular_type;
240 typedef L layout_type;
241 typedef hermitian_matrix<T, TRI, L, A> self_type;
243 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
244 using matrix_container<self_type>::operator ();
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;
255 typedef hermitian_matrix_element<self_type> reference;
257 typedef A array_type;
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;
266 // Construction and destruction
269 matrix_container<self_type> (),
270 size_ (0), data_ (0) {}
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)) {
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)) {
282 hermitian_matrix (size_type size, const array_type &data):
283 matrix_container<self_type> (),
284 size_ (size), data_ (data) {}
286 hermitian_matrix (const hermitian_matrix &m):
287 matrix_container<self_type> (),
288 size_ (m.size_), data_ (m.data_) {}
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);
300 size_type size1 () const {
304 size_type size2 () const {
310 const array_type &data () const {
314 array_type &data () {
320 void resize (size_type size, bool preserve = true) {
322 self_type temporary (size, size);
323 detail::matrix_resize_preserve<layout_type> (*this, temporary);
326 data ().resize (triangular_type::packed_size (layout_type (), size, size));
331 void resize (size_type size1, size_type size2, bool preserve = true) {
332 resize (BOOST_UBLAS_SAME (size1, size2), preserve);
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 ());
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 ());
346 // return type_traits<value_type>::real (data () [triangular_type::element (layout_type (), i, size_, i, size_)]);
348 if (triangular_type::other (i, j))
349 return data () [triangular_type::element (layout_type (), i, size_, j, size_)];
351 return type_traits<value_type>::conj (data () [triangular_type::element (layout_type (), j, size_, i, size_)]);
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_)];
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);
366 external_logic ().raise ();
367 // arbitary return value
368 return data () [triangular_type::element (layout_type (), j, size_, i, size_)];
371 if (triangular_type::other (i, j))
372 return reference (*this, i, j, data () [triangular_type::element (layout_type (), i, size_, j, size_)]);
374 return reference (*this, i, j, type_traits<value_type>::conj (data () [triangular_type::element (layout_type (), j, size_, i, size_)]));
378 // Element assignemnt
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);
386 return (data () [triangular_type::element (layout_type (), j, size_, i, size_)] = type_traits<value_type>::conj (t));
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*/();
399 std::fill (data ().begin (), data ().end (), value_type/*zero*/());
404 hermitian_matrix &operator = (const hermitian_matrix &m) {
410 hermitian_matrix &assign_temporary (hermitian_matrix &m) {
416 hermitian_matrix &operator = (const matrix_expression<AE> &ae) {
417 self_type temporary (ae);
418 return assign_temporary (temporary);
422 hermitian_matrix &assign (const matrix_expression<AE> &ae) {
423 matrix_assign<scalar_assign> (*this, ae);
428 hermitian_matrix& operator += (const matrix_expression<AE> &ae) {
429 self_type temporary (*this + ae);
430 return assign_temporary (temporary);
434 hermitian_matrix &plus_assign (const matrix_expression<AE> &ae) {
435 matrix_assign<scalar_plus_assign> (*this, ae);
440 hermitian_matrix& operator -= (const matrix_expression<AE> &ae) {
441 self_type temporary (*this - ae);
442 return assign_temporary (temporary);
446 hermitian_matrix &minus_assign (const matrix_expression<AE> &ae) {
447 matrix_assign<scalar_minus_assign> (*this, ae);
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);
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);
473 void swap (hermitian_matrix &m) {
475 std::swap (size_, m.size_);
476 data ().swap (m.data ());
480 friend void swap (hermitian_matrix &m1, hermitian_matrix &m2) {
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;
491 class const_iterator1;
493 class const_iterator2;
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;
503 const_iterator1 find1 (int /* rank */, size_type i, size_type j) const {
504 return const_iterator1 (*this, i, j);
507 iterator1 find1 (int rank, size_type i, size_type j) {
509 i = triangular_type::mutable_restrict1 (i, j);
510 return iterator1 (*this, i, j);
513 const_iterator2 find2 (int /* rank */, size_type i, size_type j) const {
514 return const_iterator2 (*this, i, j);
517 iterator2 find2 (int rank, size_type i, size_type j) {
519 j = triangular_type::mutable_restrict2 (i, j);
520 return iterator2 (*this, i, j);
523 // Iterators simply are indices.
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> {
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;
536 typedef const_iterator2 dual_iterator_type;
537 typedef const_reverse_iterator2 dual_reverse_iterator_type;
539 // Construction and destruction
542 container_const_reference<self_type> (), it1_ (), it2_ () {}
544 const_iterator1 (const self_type &m, size_type it1, size_type it2):
545 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
547 const_iterator1 (const iterator1 &it):
548 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
552 const_iterator1 &operator ++ () {
557 const_iterator1 &operator -- () {
562 const_iterator1 &operator += (difference_type n) {
567 const_iterator1 &operator -= (difference_type n) {
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_;
580 const_reference operator * () const {
581 return (*this) () (it1_, it2_);
584 const_reference operator [] (difference_type n) const {
588 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
590 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
593 const_iterator2 begin () const {
594 return (*this) ().find2 (1, it1_, 0);
597 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
600 const_iterator2 end () const {
601 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
604 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
607 const_reverse_iterator2 rbegin () const {
608 return const_reverse_iterator2 (end ());
611 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
614 const_reverse_iterator2 rend () const {
615 return const_reverse_iterator2 (begin ());
621 size_type index1 () const {
625 size_type index2 () const {
631 const_iterator1 &operator = (const const_iterator1 &it) {
632 container_const_reference<self_type>::assign (&it ());
640 bool operator == (const const_iterator1 &it) const {
641 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
642 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
643 return it1_ == it.it1_;
646 bool operator < (const const_iterator1 &it) const {
647 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
648 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
649 return it1_ < it.it1_;
659 const_iterator1 begin1 () const {
660 return find1 (0, 0, 0);
663 const_iterator1 end1 () const {
664 return find1 (0, size_, 0);
667 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
669 public container_reference<hermitian_matrix>,
670 public random_access_iterator_base<packed_random_access_iterator_tag,
671 iterator1, value_type> {
673 typedef typename hermitian_matrix::value_type value_type;
674 typedef typename hermitian_matrix::difference_type difference_type;
675 typedef typename hermitian_matrix::true_reference reference;
676 typedef typename hermitian_matrix::pointer pointer;
678 typedef iterator2 dual_iterator_type;
679 typedef reverse_iterator2 dual_reverse_iterator_type;
681 // Construction and destruction
684 container_reference<self_type> (), it1_ (), it2_ () {}
686 iterator1 (self_type &m, size_type it1, size_type it2):
687 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
691 iterator1 &operator ++ () {
696 iterator1 &operator -- () {
701 iterator1 &operator += (difference_type n) {
706 iterator1 &operator -= (difference_type n) {
711 difference_type operator - (const iterator1 &it) const {
712 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
713 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
714 return it1_ - it.it1_;
719 reference operator * () const {
720 return (*this) ().at_element (it1_, it2_);
723 reference operator [] (difference_type n) const {
727 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
729 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
732 iterator2 begin () const {
733 return (*this) ().find2 (1, it1_, 0);
736 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
739 iterator2 end () const {
740 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
743 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
746 reverse_iterator2 rbegin () const {
747 return reverse_iterator2 (end ());
750 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
753 reverse_iterator2 rend () const {
754 return reverse_iterator2 (begin ());
760 size_type index1 () const {
764 size_type index2 () const {
770 iterator1 &operator = (const iterator1 &it) {
771 container_reference<self_type>::assign (&it ());
779 bool operator == (const iterator1 &it) const {
780 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
781 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
782 return it1_ == it.it1_;
785 bool operator < (const iterator1 &it) const {
786 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
787 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
788 return it1_ < it.it1_;
795 friend class const_iterator1;
800 iterator1 begin1 () {
801 return find1 (0, 0, 0);
805 return find1 (0, size_, 0);
808 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
809 class const_iterator2:
810 public container_const_reference<hermitian_matrix>,
811 public random_access_iterator_base<packed_random_access_iterator_tag,
812 const_iterator2, value_type> {
814 typedef typename hermitian_matrix::value_type value_type;
815 typedef typename hermitian_matrix::difference_type difference_type;
816 typedef typename hermitian_matrix::const_reference reference;
817 typedef const typename hermitian_matrix::pointer pointer;
819 typedef const_iterator1 dual_iterator_type;
820 typedef const_reverse_iterator1 dual_reverse_iterator_type;
822 // Construction and destruction
825 container_const_reference<self_type> (), it1_ (), it2_ () {}
827 const_iterator2 (const self_type &m, size_type it1, size_type it2):
828 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
830 const_iterator2 (const iterator2 &it):
831 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
835 const_iterator2 &operator ++ () {
840 const_iterator2 &operator -- () {
845 const_iterator2 &operator += (difference_type n) {
850 const_iterator2 &operator -= (difference_type n) {
855 difference_type operator - (const const_iterator2 &it) const {
856 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
857 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
858 return it2_ - it.it2_;
863 const_reference operator * () const {
864 return (*this) () (it1_, it2_);
867 const_reference operator [] (difference_type n) const {
871 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
873 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
876 const_iterator1 begin () const {
877 return (*this) ().find1 (1, 0, it2_);
880 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
883 const_iterator1 end () const {
884 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
887 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
890 const_reverse_iterator1 rbegin () const {
891 return const_reverse_iterator1 (end ());
894 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
897 const_reverse_iterator1 rend () const {
898 return const_reverse_iterator1 (begin ());
904 size_type index1 () const {
908 size_type index2 () const {
914 const_iterator2 &operator = (const const_iterator2 &it) {
915 container_const_reference<self_type>::assign (&it ());
923 bool operator == (const const_iterator2 &it) const {
924 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
925 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
926 return it2_ == it.it2_;
929 bool operator < (const const_iterator2 &it) const {
930 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
931 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
932 return it2_ < it.it2_;
942 const_iterator2 begin2 () const {
943 return find2 (0, 0, 0);
946 const_iterator2 end2 () const {
947 return find2 (0, 0, size_);
950 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
952 public container_reference<hermitian_matrix>,
953 public random_access_iterator_base<packed_random_access_iterator_tag,
954 iterator2, value_type> {
956 typedef typename hermitian_matrix::value_type value_type;
957 typedef typename hermitian_matrix::difference_type difference_type;
958 typedef typename hermitian_matrix::true_reference reference;
959 typedef typename hermitian_matrix::pointer pointer;
961 typedef iterator1 dual_iterator_type;
962 typedef reverse_iterator1 dual_reverse_iterator_type;
964 // Construction and destruction
967 container_reference<self_type> (), it1_ (), it2_ () {}
969 iterator2 (self_type &m, size_type it1, size_type it2):
970 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
974 iterator2 &operator ++ () {
979 iterator2 &operator -- () {
984 iterator2 &operator += (difference_type n) {
989 iterator2 &operator -= (difference_type n) {
994 difference_type operator - (const iterator2 &it) const {
995 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
996 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
997 return it2_ - it.it2_;
1002 reference operator * () const {
1003 return (*this) ().at_element (it1_, it2_);
1006 reference operator [] (difference_type n) const {
1007 return *(*this + n);
1010 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1012 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1013 typename self_type::
1015 iterator1 begin () const {
1016 return (*this) ().find1 (1, 0, it2_);
1019 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1020 typename self_type::
1022 iterator1 end () const {
1023 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
1026 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1027 typename self_type::
1029 reverse_iterator1 rbegin () const {
1030 return reverse_iterator1 (end ());
1033 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1034 typename self_type::
1036 reverse_iterator1 rend () const {
1037 return reverse_iterator1 (begin ());
1043 size_type index1 () const {
1047 size_type index2 () const {
1053 iterator2 &operator = (const iterator2 &it) {
1054 container_reference<self_type>::assign (&it ());
1062 bool operator == (const iterator2 &it) const {
1063 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1064 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1065 return it2_ == it.it2_;
1068 bool operator < (const iterator2 &it) const {
1069 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1070 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1071 return it2_ < it.it2_;
1078 friend class const_iterator2;
1083 iterator2 begin2 () {
1084 return find2 (0, 0, 0);
1088 return find2 (0, 0, size_);
1091 // Reverse iterators
1094 const_reverse_iterator1 rbegin1 () const {
1095 return const_reverse_iterator1 (end1 ());
1098 const_reverse_iterator1 rend1 () const {
1099 return const_reverse_iterator1 (begin1 ());
1103 reverse_iterator1 rbegin1 () {
1104 return reverse_iterator1 (end1 ());
1107 reverse_iterator1 rend1 () {
1108 return reverse_iterator1 (begin1 ());
1112 const_reverse_iterator2 rbegin2 () const {
1113 return const_reverse_iterator2 (end2 ());
1116 const_reverse_iterator2 rend2 () const {
1117 return const_reverse_iterator2 (begin2 ());
1121 reverse_iterator2 rbegin2 () {
1122 return reverse_iterator2 (end2 ());
1125 reverse_iterator2 rend2 () {
1126 return reverse_iterator2 (begin2 ());
1135 // Hermitian matrix adaptor class
1136 template<class M, class TRI>
1137 class hermitian_adaptor:
1138 public matrix_expression<hermitian_adaptor<M, TRI> > {
1140 typedef hermitian_adaptor<M, TRI> self_type;
1141 typedef typename M::value_type &true_reference;
1143 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1144 using matrix_expression<self_type>::operator ();
1146 typedef const M const_matrix_type;
1147 typedef M matrix_type;
1148 typedef TRI triangular_type;
1149 typedef typename M::size_type size_type;
1150 typedef typename M::difference_type difference_type;
1151 typedef typename M::value_type value_type;
1152 typedef typename M::value_type const_reference;
1153 #ifndef BOOST_UBLAS_STRICT_HERMITIAN
1154 typedef typename boost::mpl::if_<boost::is_const<M>,
1155 typename M::value_type,
1156 typename M::reference>::type reference;
1158 typedef typename boost::mpl::if_<boost::is_const<M>,
1159 typename M::value_type,
1160 hermitian_matrix_element<self_type> >::type reference;
1162 typedef typename boost::mpl::if_<boost::is_const<M>,
1163 typename M::const_closure_type,
1164 typename M::closure_type>::type matrix_closure_type;
1165 typedef const self_type const_closure_type;
1166 typedef self_type closure_type;
1167 // Replaced by _temporary_traits to avoid type requirements on M
1168 //typedef typename M::vector_temporary_type vector_temporary_type;
1169 //typedef typename M::matrix_temporary_type matrix_temporary_type;
1170 typedef typename storage_restrict_traits<typename M::storage_category,
1171 packed_proxy_tag>::storage_category storage_category;
1172 typedef typename M::orientation_category orientation_category;
1174 // Construction and destruction
1176 hermitian_adaptor (matrix_type &data):
1177 matrix_expression<self_type> (),
1179 BOOST_UBLAS_CHECK (data_.size1 () == data_.size2 (), bad_size ());
1182 hermitian_adaptor (const hermitian_adaptor &m):
1183 matrix_expression<self_type> (),
1185 BOOST_UBLAS_CHECK (data_.size1 () == data_.size2 (), bad_size ());
1190 size_type size1 () const {
1191 return data_.size1 ();
1194 size_type size2 () const {
1195 return data_.size2 ();
1198 // Storage accessors
1200 const matrix_closure_type &data () const {
1204 matrix_closure_type &data () {
1209 #ifndef BOOST_UBLAS_PROXY_CONST_MEMBER
1211 const_reference operator () (size_type i, size_type j) const {
1212 BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1213 BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1215 // return type_traits<value_type>::real (data () (i, i));
1217 if (triangular_type::other (i, j))
1218 return data () (i, j);
1220 return type_traits<value_type>::conj (data () (j, i));
1223 reference operator () (size_type i, size_type j) {
1224 BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1225 BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1226 #ifndef BOOST_UBLAS_STRICT_HERMITIAN
1227 if (triangular_type::other (i, j))
1228 return data () (i, j);
1230 external_logic ().raise ();
1231 return conj_ = type_traits<value_type>::conj (data () (j, i));
1234 if (triangular_type::other (i, j))
1235 return reference (*this, i, j, data () (i, j));
1237 return reference (*this, i, j, type_traits<value_type>::conj (data () (j, i)));
1241 true_reference insert_element (size_type i, size_type j, value_type t) {
1242 BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1243 BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1245 // data () (i, i) = type_traits<value_type>::real (t);
1247 if (triangular_type::other (i, j))
1248 return data () (i, j) = t;
1250 return data () (j, i) = type_traits<value_type>::conj (t);
1254 reference operator () (size_type i, size_type j) {
1255 BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1256 BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1257 #ifndef BOOST_UBLAS_STRICT_HERMITIAN
1258 if (triangular_type::other (i, j))
1259 return data () (i, j);
1261 external_logic ().raise ();
1262 return conj_ = type_traits<value_type>::conj (data () (j, i));
1265 if (triangular_type::other (i, j))
1266 return reference (*this, i, j, data () (i, j));
1268 return reference (*this, i, j, type_traits<value_type>::conj (data () (j, i)));
1272 true_reference insert_element (size_type i, size_type j, value_type t) {
1273 BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1274 BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1276 // data () (i, i) = type_traits<value_type>::real (t);
1278 if (triangular_type::other (i, j))
1279 return data () (i, j) = t;
1281 return data () (j, i) = type_traits<value_type>::conj (t);
1287 hermitian_adaptor &operator = (const hermitian_adaptor &m) {
1288 matrix_assign<scalar_assign, triangular_type> (*this, m);
1292 hermitian_adaptor &assign_temporary (hermitian_adaptor &m) {
1298 hermitian_adaptor &operator = (const matrix_expression<AE> &ae) {
1299 matrix_assign<scalar_assign, triangular_type> (*this, matrix<value_type> (ae));
1304 hermitian_adaptor &assign (const matrix_expression<AE> &ae) {
1305 matrix_assign<scalar_assign, triangular_type> (*this, ae);
1310 hermitian_adaptor& operator += (const matrix_expression<AE> &ae) {
1311 matrix_assign<scalar_assign, triangular_type> (*this, matrix<value_type> (*this + ae));
1316 hermitian_adaptor &plus_assign (const matrix_expression<AE> &ae) {
1317 matrix_assign<scalar_plus_assign, triangular_type> (*this, ae);
1322 hermitian_adaptor& operator -= (const matrix_expression<AE> &ae) {
1323 matrix_assign<scalar_assign, triangular_type> (*this, matrix<value_type> (*this - ae));
1328 hermitian_adaptor &minus_assign (const matrix_expression<AE> &ae) {
1329 matrix_assign<scalar_minus_assign, triangular_type> (*this, ae);
1334 hermitian_adaptor& operator *= (const AT &at) {
1335 // Multiplication is only allowed for real scalars,
1336 // otherwise the resulting matrix isn't hermitian.
1337 // Thanks to Peter Schmitteckert for spotting this.
1338 BOOST_UBLAS_CHECK (type_traits<value_type>::imag (at) == 0, non_real ());
1339 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1344 hermitian_adaptor& operator /= (const AT &at) {
1345 // Multiplication is only allowed for real scalars,
1346 // otherwise the resulting matrix isn't hermitian.
1347 // Thanks to Peter Schmitteckert for spotting this.
1348 BOOST_UBLAS_CHECK (type_traits<value_type>::imag (at) == 0, non_real ());
1349 matrix_assign_scalar<scalar_divides_assign> (*this, at);
1353 // Closure comparison
1355 bool same_closure (const hermitian_adaptor &ha) const {
1356 return (*this).data ().same_closure (ha.data ());
1361 void swap (hermitian_adaptor &m) {
1363 matrix_swap<scalar_swap, triangular_type> (*this, m);
1366 friend void swap (hermitian_adaptor &m1, hermitian_adaptor &m2) {
1372 // Use matrix iterator
1373 typedef typename M::const_iterator1 const_subiterator1_type;
1374 typedef typename boost::mpl::if_<boost::is_const<M>,
1375 typename M::const_iterator1,
1376 typename M::iterator1>::type subiterator1_type;
1377 typedef typename M::const_iterator2 const_subiterator2_type;
1378 typedef typename boost::mpl::if_<boost::is_const<M>,
1379 typename M::const_iterator2,
1380 typename M::iterator2>::type subiterator2_type;
1383 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1384 typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
1385 typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
1386 typedef indexed_const_iterator1<self_type, dense_random_access_iterator_tag> const_iterator1;
1387 typedef indexed_const_iterator2<self_type, dense_random_access_iterator_tag> const_iterator2;
1389 class const_iterator1;
1391 class const_iterator2;
1394 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1395 typedef reverse_iterator_base1<iterator1> reverse_iterator1;
1396 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1397 typedef reverse_iterator_base2<iterator2> reverse_iterator2;
1401 const_iterator1 find1 (int rank, size_type i, size_type j) const {
1402 if (triangular_type::other (i, j)) {
1403 if (triangular_type::other (size1 (), j)) {
1404 return const_iterator1 (*this, 0, 0,
1405 data ().find1 (rank, i, j), data ().find1 (rank, size1 (), j),
1406 data ().find2 (rank, size2 (), size1 ()), data ().find2 (rank, size2 (), size1 ()));
1408 return const_iterator1 (*this, 0, 1,
1409 data ().find1 (rank, i, j), data ().find1 (rank, j, j),
1410 data ().find2 (rank, j, j), data ().find2 (rank, j, size1 ()));
1413 if (triangular_type::other (size1 (), j)) {
1414 return const_iterator1 (*this, 1, 0,
1415 data ().find1 (rank, j, j), data ().find1 (rank, size1 (), j),
1416 data ().find2 (rank, j, i), data ().find2 (rank, j, j));
1418 return const_iterator1 (*this, 1, 1,
1419 data ().find1 (rank, size1 (), size2 ()), data ().find1 (rank, size1 (), size2 ()),
1420 data ().find2 (rank, j, i), data ().find2 (rank, j, size1 ()));
1425 iterator1 find1 (int rank, size_type i, size_type j) {
1427 i = triangular_type::mutable_restrict1 (i, j);
1428 return iterator1 (*this, data ().find1 (rank, i, j));
1431 const_iterator2 find2 (int rank, size_type i, size_type j) const {
1432 if (triangular_type::other (i, j)) {
1433 if (triangular_type::other (i, size2 ())) {
1434 return const_iterator2 (*this, 1, 1,
1435 data ().find1 (rank, size2 (), size1 ()), data ().find1 (rank, size2 (), size1 ()),
1436 data ().find2 (rank, i, j), data ().find2 (rank, i, size2 ()));
1438 return const_iterator2 (*this, 1, 0,
1439 data ().find1 (rank, i, i), data ().find1 (rank, size2 (), i),
1440 data ().find2 (rank, i, j), data ().find2 (rank, i, i));
1443 if (triangular_type::other (i, size2 ())) {
1444 return const_iterator2 (*this, 0, 1,
1445 data ().find1 (rank, j, i), data ().find1 (rank, i, i),
1446 data ().find2 (rank, i, i), data ().find2 (rank, i, size2 ()));
1448 return const_iterator2 (*this, 0, 0,
1449 data ().find1 (rank, j, i), data ().find1 (rank, size2 (), i),
1450 data ().find2 (rank, size1 (), size2 ()), data ().find2 (rank, size2 (), size2 ()));
1455 iterator2 find2 (int rank, size_type i, size_type j) {
1457 j = triangular_type::mutable_restrict2 (i, j);
1458 return iterator2 (*this, data ().find2 (rank, i, j));
1461 // Iterators simply are indices.
1463 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1464 class const_iterator1:
1465 public container_const_reference<hermitian_adaptor>,
1466 public random_access_iterator_base<typename iterator_restrict_traits<
1467 typename const_subiterator1_type::iterator_category, dense_random_access_iterator_tag>::iterator_category,
1468 const_iterator1, value_type> {
1470 typedef typename const_subiterator1_type::value_type value_type;
1471 typedef typename const_subiterator1_type::difference_type difference_type;
1472 // FIXME no better way to not return the address of a temporary?
1473 // typedef typename const_subiterator1_type::reference reference;
1474 typedef typename const_subiterator1_type::value_type reference;
1475 typedef typename const_subiterator1_type::pointer pointer;
1477 typedef const_iterator2 dual_iterator_type;
1478 typedef const_reverse_iterator2 dual_reverse_iterator_type;
1480 // Construction and destruction
1483 container_const_reference<self_type> (),
1484 begin_ (-1), end_ (-1), current_ (-1),
1485 it1_begin_ (), it1_end_ (), it1_ (),
1486 it2_begin_ (), it2_end_ (), it2_ () {}
1488 const_iterator1 (const self_type &m, int begin, int end,
1489 const const_subiterator1_type &it1_begin, const const_subiterator1_type &it1_end,
1490 const const_subiterator2_type &it2_begin, const const_subiterator2_type &it2_end):
1491 container_const_reference<self_type> (m),
1492 begin_ (begin), end_ (end), current_ (begin),
1493 it1_begin_ (it1_begin), it1_end_ (it1_end), it1_ (it1_begin_),
1494 it2_begin_ (it2_begin), it2_end_ (it2_end), it2_ (it2_begin_) {
1495 if (current_ == 0 && it1_ == it1_end_)
1497 if (current_ == 1 && it2_ == it2_end_)
1499 if ((current_ == 0 && it1_ == it1_end_) ||
1500 (current_ == 1 && it2_ == it2_end_))
1502 BOOST_UBLAS_CHECK (current_ == end_ ||
1503 (current_ == 0 && it1_ != it1_end_) ||
1504 (current_ == 1 && it2_ != it2_end_), internal_logic ());
1506 // FIXME cannot compile
1507 // iterator1 does not have these members!
1509 const_iterator1 (const iterator1 &it):
1510 container_const_reference<self_type> (it ()),
1511 begin_ (it.begin_), end_ (it.end_), current_ (it.current_),
1512 it1_begin_ (it.it1_begin_), it1_end_ (it.it1_end_), it1_ (it.it1_),
1513 it2_begin_ (it.it2_begin_), it2_end_ (it.it2_end_), it2_ (it.it2_) {
1514 BOOST_UBLAS_CHECK (current_ == end_ ||
1515 (current_ == 0 && it1_ != it1_end_) ||
1516 (current_ == 1 && it2_ != it2_end_), internal_logic ());
1521 const_iterator1 &operator ++ () {
1522 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1523 if (current_ == 0) {
1524 BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1526 if (it1_ == it1_end_ && end_ == 1) {
1530 } else /* if (current_ == 1) */ {
1531 BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1533 if (it2_ == it2_end_ && end_ == 0) {
1541 const_iterator1 &operator -- () {
1542 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1543 if (current_ == 0) {
1544 if (it1_ == it1_begin_ && begin_ == 1) {
1546 BOOST_UBLAS_CHECK (it2_ != it2_begin_, internal_logic ());
1552 } else /* if (current_ == 1) */ {
1553 if (it2_ == it2_begin_ && begin_ == 0) {
1555 BOOST_UBLAS_CHECK (it1_ != it1_begin_, internal_logic ());
1565 const_iterator1 &operator += (difference_type n) {
1566 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1567 if (current_ == 0) {
1568 size_type d = (std::min) (n, it1_end_ - it1_);
1571 if (n > 0 || (end_ == 1 && it1_ == it1_end_)) {
1572 BOOST_UBLAS_CHECK (end_ == 1, external_logic ());
1573 d = (std::min) (n, it2_end_ - it2_begin_);
1574 it2_ = it2_begin_ + d;
1578 } else /* if (current_ == 1) */ {
1579 size_type d = (std::min) (n, it2_end_ - it2_);
1582 if (n > 0 || (end_ == 0 && it2_ == it2_end_)) {
1583 BOOST_UBLAS_CHECK (end_ == 0, external_logic ());
1584 d = (std::min) (n, it1_end_ - it1_begin_);
1585 it1_ = it1_begin_ + d;
1590 BOOST_UBLAS_CHECK (n == 0, external_logic ());
1594 const_iterator1 &operator -= (difference_type n) {
1595 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1596 if (current_ == 0) {
1597 size_type d = (std::min) (n, it1_ - it1_begin_);
1601 BOOST_UBLAS_CHECK (end_ == 1, external_logic ());
1602 d = (std::min) (n, it2_end_ - it2_begin_);
1603 it2_ = it2_end_ - d;
1607 } else /* if (current_ == 1) */ {
1608 size_type d = (std::min) (n, it2_ - it2_begin_);
1612 BOOST_UBLAS_CHECK (end_ == 0, external_logic ());
1613 d = (std::min) (n, it1_end_ - it1_begin_);
1614 it1_ = it1_end_ - d;
1619 BOOST_UBLAS_CHECK (n == 0, external_logic ());
1623 difference_type operator - (const const_iterator1 &it) const {
1624 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1625 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1626 BOOST_UBLAS_CHECK (it.current_ == 0 || it.current_ == 1, internal_logic ());
1627 BOOST_UBLAS_CHECK (/* begin_ == it.begin_ && */ end_ == it.end_, internal_logic ());
1628 if (current_ == 0 && it.current_ == 0) {
1629 return it1_ - it.it1_;
1630 } else if (current_ == 0 && it.current_ == 1) {
1631 if (end_ == 1 && it.end_ == 1) {
1632 return (it1_ - it.it1_end_) + (it.it2_begin_ - it.it2_);
1633 } else /* if (end_ == 0 && it.end_ == 0) */ {
1634 return (it1_ - it.it1_begin_) + (it.it2_end_ - it.it2_);
1637 } else if (current_ == 1 && it.current_ == 0) {
1638 if (end_ == 1 && it.end_ == 1) {
1639 return (it2_ - it.it2_begin_) + (it.it1_end_ - it.it1_);
1640 } else /* if (end_ == 0 && it.end_ == 0) */ {
1641 return (it2_ - it.it2_end_) + (it.it1_begin_ - it.it1_);
1643 } else /* if (current_ == 1 && it.current_ == 1) */ {
1644 return it2_ - it.it2_;
1650 const_reference operator * () const {
1651 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1652 if (current_ == 0) {
1653 BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1654 if (triangular_type::other (index1 (), index2 ()))
1657 return type_traits<value_type>::conj (*it1_);
1658 } else /* if (current_ == 1) */ {
1659 BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1660 if (triangular_type::other (index1 (), index2 ()))
1663 return type_traits<value_type>::conj (*it2_);
1667 const_reference operator [] (difference_type n) const {
1668 return *(*this + n);
1671 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1673 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1674 typename self_type::
1676 const_iterator2 begin () const {
1677 return (*this) ().find2 (1, index1 (), 0);
1680 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1681 typename self_type::
1683 const_iterator2 end () const {
1684 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1687 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1688 typename self_type::
1690 const_reverse_iterator2 rbegin () const {
1691 return const_reverse_iterator2 (end ());
1694 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1695 typename self_type::
1697 const_reverse_iterator2 rend () const {
1698 return const_reverse_iterator2 (begin ());
1704 size_type index1 () const {
1705 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1706 if (current_ == 0) {
1707 BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1708 return it1_.index1 ();
1709 } else /* if (current_ == 1) */ {
1710 BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1711 return it2_.index2 ();
1715 size_type index2 () const {
1716 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1717 if (current_ == 0) {
1718 BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1719 return it1_.index2 ();
1720 } else /* if (current_ == 1) */ {
1721 BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1722 return it2_.index1 ();
1728 const_iterator1 &operator = (const const_iterator1 &it) {
1729 container_const_reference<self_type>::assign (&it ());
1732 current_ = it.current_;
1733 it1_begin_ = it.it1_begin_;
1734 it1_end_ = it.it1_end_;
1736 it2_begin_ = it.it2_begin_;
1737 it2_end_ = it.it2_end_;
1744 bool operator == (const const_iterator1 &it) const {
1745 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1746 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1747 BOOST_UBLAS_CHECK (it.current_ == 0 || it.current_ == 1, internal_logic ());
1748 BOOST_UBLAS_CHECK (/* begin_ == it.begin_ && */ end_ == it.end_, internal_logic ());
1749 return (current_ == 0 && it.current_ == 0 && it1_ == it.it1_) ||
1750 (current_ == 1 && it.current_ == 1 && it2_ == it.it2_);
1753 bool operator < (const const_iterator1 &it) const {
1754 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1755 return it - *this > 0;
1762 const_subiterator1_type it1_begin_;
1763 const_subiterator1_type it1_end_;
1764 const_subiterator1_type it1_;
1765 const_subiterator2_type it2_begin_;
1766 const_subiterator2_type it2_end_;
1767 const_subiterator2_type it2_;
1772 const_iterator1 begin1 () const {
1773 return find1 (0, 0, 0);
1776 const_iterator1 end1 () const {
1777 return find1 (0, size1 (), 0);
1780 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1782 public container_reference<hermitian_adaptor>,
1783 public random_access_iterator_base<typename iterator_restrict_traits<
1784 typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1785 iterator1, value_type> {
1787 typedef typename subiterator1_type::value_type value_type;
1788 typedef typename subiterator1_type::difference_type difference_type;
1789 typedef typename subiterator1_type::reference reference;
1790 typedef typename subiterator1_type::pointer pointer;
1792 typedef iterator2 dual_iterator_type;
1793 typedef reverse_iterator2 dual_reverse_iterator_type;
1795 // Construction and destruction
1798 container_reference<self_type> (), it1_ () {}
1800 iterator1 (self_type &m, const subiterator1_type &it1):
1801 container_reference<self_type> (m), it1_ (it1) {}
1805 iterator1 &operator ++ () {
1810 iterator1 &operator -- () {
1815 iterator1 &operator += (difference_type n) {
1820 iterator1 &operator -= (difference_type n) {
1825 difference_type operator - (const iterator1 &it) const {
1826 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1827 return it1_ - it.it1_;
1832 reference operator * () const {
1836 reference operator [] (difference_type n) const {
1837 return *(*this + n);
1840 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1842 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1843 typename self_type::
1845 iterator2 begin () const {
1846 return (*this) ().find2 (1, index1 (), 0);
1849 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1850 typename self_type::
1852 iterator2 end () const {
1853 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1856 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1857 typename self_type::
1859 reverse_iterator2 rbegin () const {
1860 return reverse_iterator2 (end ());
1863 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1864 typename self_type::
1866 reverse_iterator2 rend () const {
1867 return reverse_iterator2 (begin ());
1873 size_type index1 () const {
1874 return it1_.index1 ();
1877 size_type index2 () const {
1878 return it1_.index2 ();
1883 iterator1 &operator = (const iterator1 &it) {
1884 container_reference<self_type>::assign (&it ());
1891 bool operator == (const iterator1 &it) const {
1892 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1893 return it1_ == it.it1_;
1896 bool operator < (const iterator1 &it) const {
1897 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1898 return it1_ < it.it1_;
1902 subiterator1_type it1_;
1904 friend class const_iterator1;
1909 iterator1 begin1 () {
1910 return find1 (0, 0, 0);
1914 return find1 (0, size1 (), 0);
1917 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1918 class const_iterator2:
1919 public container_const_reference<hermitian_adaptor>,
1920 public random_access_iterator_base<typename iterator_restrict_traits<
1921 typename const_subiterator2_type::iterator_category, dense_random_access_iterator_tag>::iterator_category,
1922 const_iterator2, value_type> {
1924 typedef typename const_subiterator2_type::value_type value_type;
1925 typedef typename const_subiterator2_type::difference_type difference_type;
1926 // FIXME no better way to not return the address of a temporary?
1927 // typedef typename const_subiterator2_type::reference reference;
1928 typedef typename const_subiterator2_type::value_type reference;
1929 typedef typename const_subiterator2_type::pointer pointer;
1931 typedef const_iterator1 dual_iterator_type;
1932 typedef const_reverse_iterator1 dual_reverse_iterator_type;
1934 // Construction and destruction
1937 container_const_reference<self_type> (),
1938 begin_ (-1), end_ (-1), current_ (-1),
1939 it1_begin_ (), it1_end_ (), it1_ (),
1940 it2_begin_ (), it2_end_ (), it2_ () {}
1942 const_iterator2 (const self_type &m, int begin, int end,
1943 const const_subiterator1_type &it1_begin, const const_subiterator1_type &it1_end,
1944 const const_subiterator2_type &it2_begin, const const_subiterator2_type &it2_end):
1945 container_const_reference<self_type> (m),
1946 begin_ (begin), end_ (end), current_ (begin),
1947 it1_begin_ (it1_begin), it1_end_ (it1_end), it1_ (it1_begin_),
1948 it2_begin_ (it2_begin), it2_end_ (it2_end), it2_ (it2_begin_) {
1949 if (current_ == 0 && it1_ == it1_end_)
1951 if (current_ == 1 && it2_ == it2_end_)
1953 if ((current_ == 0 && it1_ == it1_end_) ||
1954 (current_ == 1 && it2_ == it2_end_))
1956 BOOST_UBLAS_CHECK (current_ == end_ ||
1957 (current_ == 0 && it1_ != it1_end_) ||
1958 (current_ == 1 && it2_ != it2_end_), internal_logic ());
1960 // FIXME cannot compiler
1961 // iterator2 does not have these members!
1963 const_iterator2 (const iterator2 &it):
1964 container_const_reference<self_type> (it ()),
1965 begin_ (it.begin_), end_ (it.end_), current_ (it.current_),
1966 it1_begin_ (it.it1_begin_), it1_end_ (it.it1_end_), it1_ (it.it1_),
1967 it2_begin_ (it.it2_begin_), it2_end_ (it.it2_end_), it2_ (it.it2_) {
1968 BOOST_UBLAS_CHECK (current_ == end_ ||
1969 (current_ == 0 && it1_ != it1_end_) ||
1970 (current_ == 1 && it2_ != it2_end_), internal_logic ());
1975 const_iterator2 &operator ++ () {
1976 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1977 if (current_ == 0) {
1978 BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
1980 if (it1_ == it1_end_ && end_ == 1) {
1984 } else /* if (current_ == 1) */ {
1985 BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
1987 if (it2_ == it2_end_ && end_ == 0) {
1995 const_iterator2 &operator -- () {
1996 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
1997 if (current_ == 0) {
1998 if (it1_ == it1_begin_ && begin_ == 1) {
2000 BOOST_UBLAS_CHECK (it2_ != it2_begin_, internal_logic ());
2006 } else /* if (current_ == 1) */ {
2007 if (it2_ == it2_begin_ && begin_ == 0) {
2009 BOOST_UBLAS_CHECK (it1_ != it1_begin_, internal_logic ());
2019 const_iterator2 &operator += (difference_type n) {
2020 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2021 if (current_ == 0) {
2022 size_type d = (std::min) (n, it1_end_ - it1_);
2025 if (n > 0 || (end_ == 1 && it1_ == it1_end_)) {
2026 BOOST_UBLAS_CHECK (end_ == 1, external_logic ());
2027 d = (std::min) (n, it2_end_ - it2_begin_);
2028 it2_ = it2_begin_ + d;
2032 } else /* if (current_ == 1) */ {
2033 size_type d = (std::min) (n, it2_end_ - it2_);
2036 if (n > 0 || (end_ == 0 && it2_ == it2_end_)) {
2037 BOOST_UBLAS_CHECK (end_ == 0, external_logic ());
2038 d = (std::min) (n, it1_end_ - it1_begin_);
2039 it1_ = it1_begin_ + d;
2044 BOOST_UBLAS_CHECK (n == 0, external_logic ());
2048 const_iterator2 &operator -= (difference_type n) {
2049 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2050 if (current_ == 0) {
2051 size_type d = (std::min) (n, it1_ - it1_begin_);
2055 BOOST_UBLAS_CHECK (end_ == 1, external_logic ());
2056 d = (std::min) (n, it2_end_ - it2_begin_);
2057 it2_ = it2_end_ - d;
2061 } else /* if (current_ == 1) */ {
2062 size_type d = (std::min) (n, it2_ - it2_begin_);
2066 BOOST_UBLAS_CHECK (end_ == 0, external_logic ());
2067 d = (std::min) (n, it1_end_ - it1_begin_);
2068 it1_ = it1_end_ - d;
2073 BOOST_UBLAS_CHECK (n == 0, external_logic ());
2077 difference_type operator - (const const_iterator2 &it) const {
2078 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2079 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2080 BOOST_UBLAS_CHECK (it.current_ == 0 || it.current_ == 1, internal_logic ());
2081 BOOST_UBLAS_CHECK (/* begin_ == it.begin_ && */ end_ == it.end_, internal_logic ());
2082 if (current_ == 0 && it.current_ == 0) {
2083 return it1_ - it.it1_;
2084 } else if (current_ == 0 && it.current_ == 1) {
2085 if (end_ == 1 && it.end_ == 1) {
2086 return (it1_ - it.it1_end_) + (it.it2_begin_ - it.it2_);
2087 } else /* if (end_ == 0 && it.end_ == 0) */ {
2088 return (it1_ - it.it1_begin_) + (it.it2_end_ - it.it2_);
2091 } else if (current_ == 1 && it.current_ == 0) {
2092 if (end_ == 1 && it.end_ == 1) {
2093 return (it2_ - it.it2_begin_) + (it.it1_end_ - it.it1_);
2094 } else /* if (end_ == 0 && it.end_ == 0) */ {
2095 return (it2_ - it.it2_end_) + (it.it1_begin_ - it.it1_);
2097 } else /* if (current_ == 1 && it.current_ == 1) */ {
2098 return it2_ - it.it2_;
2104 const_reference operator * () const {
2105 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2106 if (current_ == 0) {
2107 BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
2108 if (triangular_type::other (index1 (), index2 ()))
2111 return type_traits<value_type>::conj (*it1_);
2112 } else /* if (current_ == 1) */ {
2113 BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
2114 if (triangular_type::other (index1 (), index2 ()))
2117 return type_traits<value_type>::conj (*it2_);
2121 const_reference operator [] (difference_type n) const {
2122 return *(*this + n);
2125 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2127 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2128 typename self_type::
2130 const_iterator1 begin () const {
2131 return (*this) ().find1 (1, 0, index2 ());
2134 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2135 typename self_type::
2137 const_iterator1 end () const {
2138 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2141 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2142 typename self_type::
2144 const_reverse_iterator1 rbegin () const {
2145 return const_reverse_iterator1 (end ());
2148 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2149 typename self_type::
2151 const_reverse_iterator1 rend () const {
2152 return const_reverse_iterator1 (begin ());
2158 size_type index1 () const {
2159 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2160 if (current_ == 0) {
2161 BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
2162 return it1_.index2 ();
2163 } else /* if (current_ == 1) */ {
2164 BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
2165 return it2_.index1 ();
2169 size_type index2 () const {
2170 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2171 if (current_ == 0) {
2172 BOOST_UBLAS_CHECK (it1_ != it1_end_, internal_logic ());
2173 return it1_.index1 ();
2174 } else /* if (current_ == 1) */ {
2175 BOOST_UBLAS_CHECK (it2_ != it2_end_, internal_logic ());
2176 return it2_.index2 ();
2182 const_iterator2 &operator = (const const_iterator2 &it) {
2183 container_const_reference<self_type>::assign (&it ());
2186 current_ = it.current_;
2187 it1_begin_ = it.it1_begin_;
2188 it1_end_ = it.it1_end_;
2190 it2_begin_ = it.it2_begin_;
2191 it2_end_ = it.it2_end_;
2198 bool operator == (const const_iterator2 &it) const {
2199 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2200 BOOST_UBLAS_CHECK (current_ == 0 || current_ == 1, internal_logic ());
2201 BOOST_UBLAS_CHECK (it.current_ == 0 || it.current_ == 1, internal_logic ());
2202 BOOST_UBLAS_CHECK (/* begin_ == it.begin_ && */ end_ == it.end_, internal_logic ());
2203 return (current_ == 0 && it.current_ == 0 && it1_ == it.it1_) ||
2204 (current_ == 1 && it.current_ == 1 && it2_ == it.it2_);
2207 bool operator < (const const_iterator2 &it) const {
2208 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2209 return it - *this > 0;
2216 const_subiterator1_type it1_begin_;
2217 const_subiterator1_type it1_end_;
2218 const_subiterator1_type it1_;
2219 const_subiterator2_type it2_begin_;
2220 const_subiterator2_type it2_end_;
2221 const_subiterator2_type it2_;
2226 const_iterator2 begin2 () const {
2227 return find2 (0, 0, 0);
2230 const_iterator2 end2 () const {
2231 return find2 (0, 0, size2 ());
2234 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2236 public container_reference<hermitian_adaptor>,
2237 public random_access_iterator_base<typename iterator_restrict_traits<
2238 typename subiterator2_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
2239 iterator2, value_type> {
2241 typedef typename subiterator2_type::value_type value_type;
2242 typedef typename subiterator2_type::difference_type difference_type;
2243 typedef typename subiterator2_type::reference reference;
2244 typedef typename subiterator2_type::pointer pointer;
2246 typedef iterator1 dual_iterator_type;
2247 typedef reverse_iterator1 dual_reverse_iterator_type;
2249 // Construction and destruction
2252 container_reference<self_type> (), it2_ () {}
2254 iterator2 (self_type &m, const subiterator2_type &it2):
2255 container_reference<self_type> (m), it2_ (it2) {}
2259 iterator2 &operator ++ () {
2264 iterator2 &operator -- () {
2269 iterator2 &operator += (difference_type n) {
2274 iterator2 &operator -= (difference_type n) {
2279 difference_type operator - (const iterator2 &it) const {
2280 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2281 return it2_ - it.it2_;
2286 reference operator * () const {
2290 reference operator [] (difference_type n) const {
2291 return *(*this + n);
2294 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2296 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2297 typename self_type::
2299 iterator1 begin () const {
2300 return (*this) ().find1 (1, 0, index2 ());
2303 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2304 typename self_type::
2306 iterator1 end () const {
2307 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2310 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2311 typename self_type::
2313 reverse_iterator1 rbegin () const {
2314 return reverse_iterator1 (end ());
2317 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2318 typename self_type::
2320 reverse_iterator1 rend () const {
2321 return reverse_iterator1 (begin ());
2327 size_type index1 () const {
2328 return it2_.index1 ();
2331 size_type index2 () const {
2332 return it2_.index2 ();
2337 iterator2 &operator = (const iterator2 &it) {
2338 container_reference<self_type>::assign (&it ());
2345 bool operator == (const iterator2 &it) const {
2346 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2347 return it2_ == it.it2_;
2350 bool operator < (const iterator2 &it) const {
2351 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2352 return it2_ < it.it2_;
2356 subiterator2_type it2_;
2358 friend class const_iterator2;
2363 iterator2 begin2 () {
2364 return find2 (0, 0, 0);
2368 return find2 (0, 0, size2 ());
2371 // Reverse iterators
2374 const_reverse_iterator1 rbegin1 () const {
2375 return const_reverse_iterator1 (end1 ());
2378 const_reverse_iterator1 rend1 () const {
2379 return const_reverse_iterator1 (begin1 ());
2383 reverse_iterator1 rbegin1 () {
2384 return reverse_iterator1 (end1 ());
2387 reverse_iterator1 rend1 () {
2388 return reverse_iterator1 (begin1 ());
2392 const_reverse_iterator2 rbegin2 () const {
2393 return const_reverse_iterator2 (end2 ());
2396 const_reverse_iterator2 rend2 () const {
2397 return const_reverse_iterator2 (begin2 ());
2401 reverse_iterator2 rbegin2 () {
2402 return reverse_iterator2 (end2 ());
2405 reverse_iterator2 rend2 () {
2406 return reverse_iterator2 (begin2 ());
2410 matrix_closure_type data_;
2411 static value_type conj_;
2414 template<class M, class TRI>
2415 typename hermitian_adaptor<M, TRI>::value_type hermitian_adaptor<M, TRI>::conj_;
2417 // Specialization for temporary_traits
2418 template <class M, class TRI>
2419 struct vector_temporary_traits< hermitian_adaptor<M, TRI> >
2420 : vector_temporary_traits< M > {} ;
2421 template <class M, class TRI>
2422 struct vector_temporary_traits< const hermitian_adaptor<M, TRI> >
2423 : vector_temporary_traits< M > {} ;
2425 template <class M, class TRI>
2426 struct matrix_temporary_traits< hermitian_adaptor<M, TRI> >
2427 : matrix_temporary_traits< M > {} ;
2428 template <class M, class TRI>
2429 struct matrix_temporary_traits< const hermitian_adaptor<M, TRI> >
2430 : matrix_temporary_traits< M > {} ;