2 // Copyright (c) 2000-2013
3 // Joerg Walter, Mathias Koch, Athanasios Iliopoulos
5 // Distributed under the Boost Software License, Version 1.0. (See
6 // accompanying file LICENSE_1_0.txt or copy at
7 // http://www.boost.org/LICENSE_1_0.txt)
9 // The authors gratefully acknowledge the support of
10 // GeNeSys mbH & Co. KG in producing this work.
13 #ifndef _BOOST_UBLAS_BANDED_
14 #define _BOOST_UBLAS_BANDED_
16 #include <boost/numeric/ublas/matrix.hpp>
17 #include <boost/numeric/ublas/detail/temporary.hpp>
19 // Iterators based on ideas of Jeremy Siek
21 namespace boost { namespace numeric { namespace ublas {
28 /** \brief A helper for band_matrix indexing.
30 * The indexing happens as per the netlib description: http://www.netlib.org/lapack/lug/node124.html.
31 * In the case of a row_major matrix a different approach is followed;
33 template <class LayoutType>
34 class banded_indexing { };
36 /** \brief A helper for indexing column major banded matrices.
40 class banded_indexing<column_major_tag> {
44 BOOST_UBLAS_INLINE static T size(T /*size1*/, T size2) {
49 // BOOST_UBLAS_INLINE static bool valid_index(T size1, T /*size2*/, T lower, T upper, T i, T j) {
50 // return (upper+i >= j) && i <= std::min(size1 - 1, j + lower); // upper + i is used by get_index. Maybe find a way to consolidate the operations to increase performance
54 BOOST_UBLAS_INLINE static T get_index(T /*size1*/, T size2, T lower, T upper, T i, T j) {
55 return column_major::element (upper + i - j, lower + 1 + upper, j, size2);
59 /** \brief A helper for indexing row major banded matrices.
63 class banded_indexing<row_major_tag> {
67 BOOST_UBLAS_INLINE static T size(T size1, T /*size2*/) {
72 // BOOST_UBLAS_INLINE static bool valid_index(T /*size1*/, T size2, T lower, T upper, T i, T j) {
73 // return (lower+j >= i) && j <= std::min(size2 - 1, i + upper); // lower + j is used by get_index. Maybe find a way to consolidate the operations to increase performance
77 BOOST_UBLAS_INLINE static T get_index(T size1, T /*size2*/, T lower, T upper, T i, T j) {
78 return row_major::element (i, size1, lower + j - i, lower + 1 + upper);
84 /** \brief A banded matrix of values of type \c T.
86 * For a \f$(mxn)\f$-dimensional banded matrix with \f$l\f$ lower and \f$u\f$ upper diagonals and
87 * \f$0 \leq i < m\f$ and \f$0 \leq j < n\f$, if \f$i>j+l\f$ or \f$i<j-u\f$ then \f$b_{i,j}=0\f$.
88 * The default storage for banded matrices is packed. Orientation and storage can also be specified.
89 * Default is \c row_major and and unbounded_array. It is \b not required by the storage to initialize
90 * elements of the matrix.
92 * \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
93 * \tparam L the storage organization. It can be either \c row_major or \c column_major. Default is \c row_major
94 * \tparam A the type of Storage array. Default is \c unbounded_array
96 template<class T, class L, class A>
98 public matrix_container<banded_matrix<T, L, A> > {
101 typedef L layout_type;
102 typedef banded_matrix<T, L, A> self_type;
107 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
108 using matrix_container<self_type>::operator ();
110 typedef typename A::size_type size_type;
111 typedef typename A::difference_type difference_type;
112 typedef T value_type;
113 typedef const T &const_reference;
114 typedef T &reference;
115 typedef A array_type;
116 typedef const matrix_reference<const self_type> const_closure_type;
117 typedef matrix_reference<self_type> closure_type;
118 typedef vector<T, A> vector_temporary_type;
119 typedef matrix<T, L, A> matrix_temporary_type; // general sub-matrix
120 typedef packed_tag storage_category;
121 typedef typename L::orientation_category orientation_category;
126 // Construction and destruction
129 matrix_container<self_type> (),
130 size1_ (0), size2_ (0),
131 lower_ (0), upper_ (0), data_ (0) {}
133 banded_matrix (size_type size1, size_type size2, size_type lower = 0, size_type upper = 0):
134 matrix_container<self_type> (),
135 size1_ (size1), size2_ (size2),
136 lower_ (lower), upper_ (upper),
137 #if defined(BOOST_UBLAS_OWN_BANDED) || (BOOST_UBLAS_LEGACY_BANDED)
138 data_ ((std::max) (size1, size2) * (lower + 1 + upper))
140 data_ ( hidden::banded_indexing<orientation_category>::size(size1, size2) * (lower + 1 + upper)) // This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
145 banded_matrix (size_type size1, size_type size2, size_type lower, size_type upper, const array_type &data):
146 matrix_container<self_type> (),
147 size1_ (size1), size2_ (size2),
148 lower_ (lower), upper_ (upper), data_ (data) {}
150 banded_matrix (const banded_matrix &m):
151 matrix_container<self_type> (),
152 size1_ (m.size1_), size2_ (m.size2_),
153 lower_ (m.lower_), upper_ (m.upper_), data_ (m.data_) {}
156 banded_matrix (const matrix_expression<AE> &ae, size_type lower = 0, size_type upper = 0):
157 matrix_container<self_type> (),
158 size1_ (ae ().size1 ()), size2_ (ae ().size2 ()),
159 lower_ (lower), upper_ (upper),
160 #if defined(BOOST_UBLAS_OWN_BANDED) || (BOOST_UBLAS_LEGACY_BANDED)
161 data_ ((std::max) (size1_, size2_) * (lower_ + 1 + upper_))
163 data_ ( hidden::banded_indexing<orientation_category>::size(size1_, size2_) * (lower_ + 1 + upper_)) // This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
166 matrix_assign<scalar_assign> (*this, ae);
171 size_type size1 () const {
175 size_type size2 () const {
179 size_type lower () const {
183 size_type upper () const {
189 const array_type &data () const {
193 array_type &data () {
197 #if !defined (BOOST_UBLAS_OWN_BANDED)||(BOOST_UBLAS_LEGACY_BANDED)
199 bool is_element_in_band(size_type i, size_type j) const{
200 //return (upper_+i >= j) && i <= std::min(size1() - 1, j + lower_); // We don't need to check if i is outside because it is checked anyway in the accessors.
201 return (upper_+i >= j) && i <= ( j + lower_); // Essentially this band has "infinite" positive dimensions
206 void resize (size_type size1, size_type size2, size_type lower = 0, size_type upper = 0, bool preserve = true) {
208 self_type temporary (size1, size2, lower, upper);
209 detail::matrix_resize_preserve<layout_type> (*this, temporary);
212 data ().resize ((std::max) (size1, size2) * (lower + 1 + upper));
221 void resize_packed_preserve (size_type size1, size_type size2, size_type lower = 0, size_type upper = 0) {
226 data ().resize ((std::max) (size1, size2) * (lower + 1 + upper), value_type ());
231 const_reference operator () (size_type i, size_type j) const {
232 BOOST_UBLAS_CHECK (i < size1_, bad_index ());
233 BOOST_UBLAS_CHECK (j < size2_, bad_index ());
234 #ifdef BOOST_UBLAS_OWN_BANDED
235 const size_type k = (std::max) (i, j);
236 const size_type l = lower_ + j - i;
237 if (k < (std::max) (size1_, size2_) && // TODO: probably use BOOST_UBLAS_CHECK here instead of if
238 l < lower_ + 1 + upper_)
239 return data () [layout_type::element (k, (std::max) (size1_, size2_),
240 l, lower_ + 1 + upper_)];
241 #elif BOOST_UBLAS_LEGACY_BANDED // Prior to version: TODO: add version this is actually incorporated in
242 const size_type k = j;
243 const size_type l = upper_ + i - j;
245 l < lower_ + 1 + upper_)
246 return data () [layout_type::element (k, size2_,
247 l, lower_ + 1 + upper_)];
249 // This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
250 if ( is_element_in_band( i, j) ) {
251 return data () [hidden::banded_indexing<orientation_category>::get_index(size1_, size2_, lower_, upper_, i, j)];
258 reference at_element (size_type i, size_type j) {
259 BOOST_UBLAS_CHECK (i < size1_, bad_index ());
260 BOOST_UBLAS_CHECK (j < size2_, bad_index ());
261 #ifdef BOOST_UBLAS_OWN_BANDED
262 const size_type k = (std::max) (i, j);
263 const size_type l = lower_ + j - i; // TODO: Don't we need an if or BOOST_UBLAS_CHECK HERE?
264 return data () [layout_type::element (k, (std::max) (size1_, size2_),
265 l, lower_ + 1 + upper_)];
266 #elif BOOST_UBLAS_LEGACY_BANDED // Prior to version: TODO: add version this is actually incorporated in
267 const size_type k = j;
268 const size_type l = upper_ + i - j;
270 l < lower_ + 1 + upper_) ) {
271 bad_index ().raise ();
274 return data () [layout_type::element (k, size2_,
275 l, lower_ + 1 + upper_)];
277 // This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
278 BOOST_UBLAS_CHECK(is_element_in_band( i, j) , bad_index());
279 return data () [hidden::banded_indexing<orientation_category>::get_index(size1_, size2_, lower_, upper_, i, j)];
283 reference operator () (size_type i, size_type j) {
284 BOOST_UBLAS_CHECK (i < size1_, bad_index ());
285 BOOST_UBLAS_CHECK (j < size2_, bad_index ());
286 #ifdef BOOST_UBLAS_OWN_BANDED
287 const size_type k = (std::max) (i, j);
288 const size_type l = lower_ + j - i;
289 if (! (k < (std::max) (size1_, size2_) && // TODO: probably use BOOST_UBLAS_CHECK here instead of if
290 l < lower_ + 1 + upper_) ) {
291 bad_index ().raise ();
294 return data () [layout_type::element (k, (std::max) (size1_, size2_),
295 l, lower_ + 1 + upper_)];
296 #elif BOOST_UBLAS_LEGACY_BANDED // Prior to version: TODO: add version this is actually incorporated in
297 const size_type k = j;
298 const size_type l = upper_ + i - j;
300 l < lower_ + 1 + upper_) ) {
301 bad_index ().raise ();
304 return data () [layout_type::element (k, size2_,
305 l, lower_ + 1 + upper_)];
307 // This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
308 BOOST_UBLAS_CHECK( is_element_in_band( i, j) , bad_index());
309 return data () [hidden::banded_indexing<orientation_category>::get_index(size1_, size2_, lower_, upper_, i, j)];
314 // Element assignment
316 reference insert_element (size_type i, size_type j, const_reference t) {
317 return (operator () (i, j) = t);
320 void erase_element (size_type i, size_type j) {
321 operator () (i, j) = value_type/*zero*/();
327 std::fill (data ().begin (), data ().end (), value_type/*zero*/());
332 banded_matrix &operator = (const banded_matrix &m) {
341 banded_matrix &assign_temporary (banded_matrix &m) {
347 banded_matrix &operator = (const matrix_expression<AE> &ae) {
348 self_type temporary (ae, lower_, upper_);
349 return assign_temporary (temporary);
353 banded_matrix &assign (const matrix_expression<AE> &ae) {
354 matrix_assign<scalar_assign> (*this, ae);
359 banded_matrix& operator += (const matrix_expression<AE> &ae) {
360 self_type temporary (*this + ae, lower_, upper_);
361 return assign_temporary (temporary);
365 banded_matrix &plus_assign (const matrix_expression<AE> &ae) {
366 matrix_assign<scalar_plus_assign> (*this, ae);
371 banded_matrix& operator -= (const matrix_expression<AE> &ae) {
372 self_type temporary (*this - ae, lower_, upper_);
373 return assign_temporary (temporary);
377 banded_matrix &minus_assign (const matrix_expression<AE> &ae) {
378 matrix_assign<scalar_minus_assign> (*this, ae);
383 banded_matrix& operator *= (const AT &at) {
384 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
389 banded_matrix& operator /= (const AT &at) {
390 matrix_assign_scalar<scalar_divides_assign> (*this, at);
396 void swap (banded_matrix &m) {
398 std::swap (size1_, m.size1_);
399 std::swap (size2_, m.size2_);
400 std::swap (lower_, m.lower_);
401 std::swap (upper_, m.upper_);
402 data ().swap (m.data ());
406 friend void swap (banded_matrix &m1, banded_matrix &m2) {
411 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
412 typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
413 typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
414 typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
415 typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
417 class const_iterator1;
419 class const_iterator2;
422 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
423 typedef reverse_iterator_base1<iterator1> reverse_iterator1;
424 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
425 typedef reverse_iterator_base2<iterator2> reverse_iterator2;
429 const_iterator1 find1 (int rank, size_type i, size_type j) const {
431 size_type lower_i = (std::max) (difference_type (j - upper_), difference_type (0));
432 i = (std::max) (i, lower_i);
433 size_type upper_i = (std::min) (j + 1 + lower_, size1_);
434 i = (std::min) (i, upper_i);
436 return const_iterator1 (*this, i, j);
439 iterator1 find1 (int rank, size_type i, size_type j) {
441 size_type lower_i = (std::max) (difference_type (j - upper_), difference_type (0));
442 i = (std::max) (i, lower_i);
443 size_type upper_i = (std::min) (j + 1 + lower_, size1_);
444 i = (std::min) (i, upper_i);
446 return iterator1 (*this, i, j);
449 const_iterator2 find2 (int rank, size_type i, size_type j) const {
451 size_type lower_j = (std::max) (difference_type (i - lower_), difference_type (0));
452 j = (std::max) (j, lower_j);
453 size_type upper_j = (std::min) (i + 1 + upper_, size2_);
454 j = (std::min) (j, upper_j);
456 return const_iterator2 (*this, i, j);
459 iterator2 find2 (int rank, size_type i, size_type j) {
461 size_type lower_j = (std::max) (difference_type (i - lower_), difference_type (0));
462 j = (std::max) (j, lower_j);
463 size_type upper_j = (std::min) (i + 1 + upper_, size2_);
464 j = (std::min) (j, upper_j);
466 return iterator2 (*this, i, j);
469 // Iterators simply are indices.
471 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
472 class const_iterator1:
473 public container_const_reference<banded_matrix>,
474 public random_access_iterator_base<packed_random_access_iterator_tag,
475 const_iterator1, value_type> {
477 typedef typename banded_matrix::value_type value_type;
478 typedef typename banded_matrix::difference_type difference_type;
479 typedef typename banded_matrix::const_reference reference;
480 typedef const typename banded_matrix::pointer pointer;
482 typedef const_iterator2 dual_iterator_type;
483 typedef const_reverse_iterator2 dual_reverse_iterator_type;
485 // Construction and destruction
488 container_const_reference<self_type> (), it1_ (), it2_ () {}
490 const_iterator1 (const self_type &m, size_type it1, size_type it2):
491 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
493 const_iterator1 (const iterator1 &it):
494 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
498 const_iterator1 &operator ++ () {
503 const_iterator1 &operator -- () {
508 const_iterator1 &operator += (difference_type n) {
513 const_iterator1 &operator -= (difference_type n) {
518 difference_type operator - (const const_iterator1 &it) const {
519 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
520 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
521 return it1_ - it.it1_;
526 const_reference operator * () const {
527 return (*this) () (it1_, it2_);
530 const_reference operator [] (difference_type n) const {
534 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
536 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
539 const_iterator2 begin () const {
540 return (*this) ().find2 (1, it1_, 0);
543 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
546 const_iterator2 cbegin () const {
550 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
553 const_iterator2 end () const {
554 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
557 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
560 const_iterator2 cend () const {
564 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
567 const_reverse_iterator2 rbegin () const {
568 return const_reverse_iterator2 (end ());
571 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
574 const_reverse_iterator2 crbegin () const {
578 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
581 const_reverse_iterator2 rend () const {
582 return const_reverse_iterator2 (begin ());
585 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
588 const_reverse_iterator2 crend () const {
595 size_type index1 () const {
599 size_type index2 () const {
605 const_iterator1 &operator = (const const_iterator1 &it) {
606 container_const_reference<self_type>::assign (&it ());
614 bool operator == (const const_iterator1 &it) const {
615 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
616 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
617 return it1_ == it.it1_;
620 bool operator < (const const_iterator1 &it) const {
621 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
622 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
623 return it1_ < it.it1_;
633 const_iterator1 begin1 () const {
634 return find1 (0, 0, 0);
637 const_iterator1 cbegin1 () const {
641 const_iterator1 end1 () const {
642 return find1 (0, size1_, 0);
645 const_iterator1 cend1 () const {
649 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
651 public container_reference<banded_matrix>,
652 public random_access_iterator_base<packed_random_access_iterator_tag,
653 iterator1, value_type> {
655 typedef typename banded_matrix::value_type value_type;
656 typedef typename banded_matrix::difference_type difference_type;
657 typedef typename banded_matrix::reference reference;
658 typedef typename banded_matrix::pointer pointer;
660 typedef iterator2 dual_iterator_type;
661 typedef reverse_iterator2 dual_reverse_iterator_type;
663 // Construction and destruction
666 container_reference<self_type> (), it1_ (), it2_ () {}
668 iterator1 (self_type &m, size_type it1, size_type it2):
669 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
673 iterator1 &operator ++ () {
678 iterator1 &operator -- () {
683 iterator1 &operator += (difference_type n) {
688 iterator1 &operator -= (difference_type n) {
693 difference_type operator - (const iterator1 &it) const {
694 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
695 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
696 return it1_ - it.it1_;
701 reference operator * () const {
702 return (*this) ().at_element (it1_, it2_);
705 reference operator [] (difference_type n) const {
709 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
711 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
714 iterator2 begin () const {
715 return (*this) ().find2 (1, it1_, 0);
718 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
721 iterator2 end () const {
722 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
726 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
729 reverse_iterator2 rbegin () const {
730 return reverse_iterator2 (end ());
733 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
736 reverse_iterator2 rend () const {
737 return reverse_iterator2 (begin ());
743 size_type index1 () const {
747 size_type index2 () const {
753 iterator1 &operator = (const iterator1 &it) {
754 container_reference<self_type>::assign (&it ());
762 bool operator == (const iterator1 &it) const {
763 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
764 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
765 return it1_ == it.it1_;
768 bool operator < (const iterator1 &it) const {
769 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
770 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
771 return it1_ < it.it1_;
778 friend class const_iterator1;
783 iterator1 begin1 () {
784 return find1 (0, 0, 0);
788 return find1 (0, size1_, 0);
791 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
792 class const_iterator2:
793 public container_const_reference<banded_matrix>,
794 public random_access_iterator_base<packed_random_access_iterator_tag,
795 const_iterator2, value_type> {
797 typedef typename banded_matrix::value_type value_type;
798 typedef typename banded_matrix::difference_type difference_type;
799 typedef typename banded_matrix::const_reference reference;
800 typedef const typename banded_matrix::pointer pointer;
802 typedef const_iterator1 dual_iterator_type;
803 typedef const_reverse_iterator1 dual_reverse_iterator_type;
805 // Construction and destruction
808 container_const_reference<self_type> (), it1_ (), it2_ () {}
810 const_iterator2 (const self_type &m, size_type it1, size_type it2):
811 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
813 const_iterator2 (const iterator2 &it):
814 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
818 const_iterator2 &operator ++ () {
823 const_iterator2 &operator -- () {
828 const_iterator2 &operator += (difference_type n) {
833 const_iterator2 &operator -= (difference_type n) {
838 difference_type operator - (const const_iterator2 &it) const {
839 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
840 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
841 return it2_ - it.it2_;
846 const_reference operator * () const {
847 return (*this) () (it1_, it2_);
850 const_reference operator [] (difference_type n) const {
854 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
856 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
859 const_iterator1 begin () const {
860 return (*this) ().find1 (1, 0, it2_);
863 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
866 const_iterator1 cbegin () const {
870 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
873 const_iterator1 end () const {
874 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
877 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
880 const_iterator1 cend () const {
885 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
888 const_reverse_iterator1 rbegin () const {
889 return const_reverse_iterator1 (end ());
892 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
895 const_reverse_iterator1 crbegin () const {
899 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
902 const_reverse_iterator1 rend () const {
903 return const_reverse_iterator1 (begin ());
906 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
909 const_reverse_iterator1 crend () const {
916 size_type index1 () const {
920 size_type index2 () const {
926 const_iterator2 &operator = (const const_iterator2 &it) {
927 container_const_reference<self_type>::assign (&it ());
935 bool operator == (const const_iterator2 &it) const {
936 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
937 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
938 return it2_ == it.it2_;
941 bool operator < (const const_iterator2 &it) const {
942 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
943 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
944 return it2_ < it.it2_;
954 const_iterator2 begin2 () const {
955 return find2 (0, 0, 0);
958 const_iterator2 cbegin2 () const {
962 const_iterator2 end2 () const {
963 return find2 (0, 0, size2_);
966 const_iterator2 cend2 () const {
970 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
972 public container_reference<banded_matrix>,
973 public random_access_iterator_base<packed_random_access_iterator_tag,
974 iterator2, value_type> {
976 typedef typename banded_matrix::value_type value_type;
977 typedef typename banded_matrix::difference_type difference_type;
978 typedef typename banded_matrix::reference reference;
979 typedef typename banded_matrix::pointer pointer;
981 typedef iterator1 dual_iterator_type;
982 typedef reverse_iterator1 dual_reverse_iterator_type;
984 // Construction and destruction
987 container_reference<self_type> (), it1_ (), it2_ () {}
989 iterator2 (self_type &m, size_type it1, size_type it2):
990 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
994 iterator2 &operator ++ () {
999 iterator2 &operator -- () {
1004 iterator2 &operator += (difference_type n) {
1009 iterator2 &operator -= (difference_type n) {
1014 difference_type operator - (const iterator2 &it) const {
1015 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1016 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1017 return it2_ - it.it2_;
1022 reference operator * () const {
1023 return (*this) ().at_element (it1_, it2_);
1026 reference operator [] (difference_type n) const {
1027 return *(*this + n);
1030 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1032 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1033 typename self_type::
1035 iterator1 begin () const {
1036 return (*this) ().find1 (1, 0, it2_);
1039 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1040 typename self_type::
1042 iterator1 end () const {
1043 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
1046 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1047 typename self_type::
1049 reverse_iterator1 rbegin () const {
1050 return reverse_iterator1 (end ());
1053 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1054 typename self_type::
1056 reverse_iterator1 rend () const {
1057 return reverse_iterator1 (begin ());
1063 size_type index1 () const {
1067 size_type index2 () const {
1073 iterator2 &operator = (const iterator2 &it) {
1074 container_reference<self_type>::assign (&it ());
1082 bool operator == (const iterator2 &it) const {
1083 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1084 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1085 return it2_ == it.it2_;
1088 bool operator < (const iterator2 &it) const {
1089 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1090 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1091 return it2_ < it.it2_;
1098 friend class const_iterator2;
1103 iterator2 begin2 () {
1104 return find2 (0, 0, 0);
1108 return find2 (0, 0, size2_);
1111 // Reverse iterators
1114 const_reverse_iterator1 rbegin1 () const {
1115 return const_reverse_iterator1 (end1 ());
1118 const_reverse_iterator1 crbegin1 () const {
1122 const_reverse_iterator1 rend1 () const {
1123 return const_reverse_iterator1 (begin1 ());
1126 const_reverse_iterator1 crend1 () const {
1131 reverse_iterator1 rbegin1 () {
1132 return reverse_iterator1 (end1 ());
1135 reverse_iterator1 rend1 () {
1136 return reverse_iterator1 (begin1 ());
1140 const_reverse_iterator2 rbegin2 () const {
1141 return const_reverse_iterator2 (end2 ());
1144 const_reverse_iterator2 crbegin2 () const {
1148 const_reverse_iterator2 rend2 () const {
1149 return const_reverse_iterator2 (begin2 ());
1152 const_reverse_iterator2 crend2 () const {
1157 reverse_iterator2 rbegin2 () {
1158 return reverse_iterator2 (end2 ());
1161 reverse_iterator2 rend2 () {
1162 return reverse_iterator2 (begin2 ());
1171 typedef const value_type const_value_type;
1172 static const_value_type zero_;
1175 template<class T, class L, class A>
1176 typename banded_matrix<T, L, A>::const_value_type banded_matrix<T, L, A>::zero_ = value_type/*zero*/();
1179 /** \brief A diagonal matrix of values of type \c T, which is a specialization of a banded matrix
1181 * For a \f$(m\times m)\f$-dimensional diagonal matrix, \f$0 \leq i < m\f$ and \f$0 \leq j < m\f$,
1182 * if \f$i\neq j\f$ then \f$b_{i,j}=0\f$. The default storage for diagonal matrices is packed.
1183 * Orientation and storage can also be specified. Default is \c row major \c unbounded_array.
1185 * As a specialization of a banded matrix, the constructor of the diagonal matrix creates
1186 * a banded matrix with 0 upper and lower diagonals around the main diagonal and the matrix is
1187 * obviously a square matrix. Operations are optimized based on these 2 assumptions. It is
1188 * \b not required by the storage to initialize elements of the matrix.
1190 * \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
1191 * \tparam L the storage organization. It can be either \c row_major or \c column_major. Default is \c row_major
1192 * \tparam A the type of Storage array. Default is \c unbounded_array
1194 template<class T, class L, class A>
1195 class diagonal_matrix:
1196 public banded_matrix<T, L, A> {
1198 typedef typename A::size_type size_type;
1199 typedef banded_matrix<T, L, A> matrix_type;
1200 typedef A array_type;
1202 // Construction and destruction
1207 diagonal_matrix (size_type size):
1208 matrix_type (size, size) {}
1210 diagonal_matrix (size_type size, const array_type& data):
1211 matrix_type (size, size, 0, 0, data) {}
1213 diagonal_matrix (size_type size1, size_type size2):
1214 matrix_type (size1, size2) {}
1217 diagonal_matrix (const matrix_expression<AE> &ae):
1220 ~diagonal_matrix () {}
1224 diagonal_matrix &operator = (const diagonal_matrix &m) {
1225 matrix_type::operator = (m);
1230 diagonal_matrix &operator = (const matrix_expression<AE> &ae) {
1231 matrix_type::operator = (ae);
1236 /** \brief A banded matrix adaptator: convert a any matrix into a banded matrix expression
1238 * For a \f$(m\times n)\f$-dimensional matrix, the \c banded_adaptor will provide a banded matrix
1239 * with \f$l\f$ lower and \f$u\f$ upper diagonals and \f$0 \leq i < m\f$ and \f$0 \leq j < n\f$,
1240 * if \f$i>j+l\f$ or \f$i<j-u\f$ then \f$b_{i,j}=0\f$.
1242 * Storage and location are based on those of the underlying matrix. This is important because
1243 * a \c banded_adaptor does not copy the matrix data to a new place. Therefore, modifying values
1244 * in a \c banded_adaptor matrix will also modify the underlying matrix too.
1246 * \tparam M the type of matrix used to generate a banded matrix
1249 class banded_adaptor:
1250 public matrix_expression<banded_adaptor<M> > {
1252 typedef banded_adaptor<M> self_type;
1255 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1256 using matrix_expression<self_type>::operator ();
1258 typedef const M const_matrix_type;
1259 typedef M matrix_type;
1260 typedef typename M::size_type size_type;
1261 typedef typename M::difference_type difference_type;
1262 typedef typename M::value_type value_type;
1263 typedef typename M::const_reference const_reference;
1264 typedef typename boost::mpl::if_<boost::is_const<M>,
1265 typename M::const_reference,
1266 typename M::reference>::type reference;
1267 typedef typename boost::mpl::if_<boost::is_const<M>,
1268 typename M::const_closure_type,
1269 typename M::closure_type>::type matrix_closure_type;
1270 typedef const self_type const_closure_type;
1271 typedef self_type closure_type;
1272 // Replaced by _temporary_traits to avoid type requirements on M
1273 //typedef typename M::vector_temporary_type vector_temporary_type;
1274 //typedef typename M::matrix_temporary_type matrix_temporary_type;
1275 typedef typename storage_restrict_traits<typename M::storage_category,
1276 packed_proxy_tag>::storage_category storage_category;
1277 typedef typename M::orientation_category orientation_category;
1279 // Construction and destruction
1281 banded_adaptor (matrix_type &data, size_type lower = 0, size_type upper = 0):
1282 matrix_expression<self_type> (),
1283 data_ (data), lower_ (lower), upper_ (upper) {}
1285 banded_adaptor (const banded_adaptor &m):
1286 matrix_expression<self_type> (),
1287 data_ (m.data_), lower_ (m.lower_), upper_ (m.upper_) {}
1291 size_type size1 () const {
1292 return data_.size1 ();
1295 size_type size2 () const {
1296 return data_.size2 ();
1299 size_type lower () const {
1303 size_type upper () const {
1307 // Storage accessors
1309 const matrix_closure_type &data () const {
1313 matrix_closure_type &data () {
1317 #if !defined (BOOST_UBLAS_OWN_BANDED)||(BOOST_UBLAS_LEGACY_BANDED)
1319 bool is_element_in_band(size_type i, size_type j) const{
1320 //return (upper_+i >= j) && i <= std::min(size1() - 1, j + lower_); // We don't need to check if i is outside because it is checked anyway in the accessors.
1321 return (upper_+i >= j) && i <= ( j + lower_); // Essentially this band has "infinite" positive dimensions
1326 #ifndef BOOST_UBLAS_PROXY_CONST_MEMBER
1328 const_reference operator () (size_type i, size_type j) const {
1329 BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1330 BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1331 #ifdef BOOST_UBLAS_OWN_BANDED
1332 size_type k = (std::max) (i, j);
1333 size_type l = lower_ + j - i;
1334 if (k < (std::max) (size1 (), size2 ()) &&
1335 l < lower_ + 1 + upper_)
1336 return data () (i, j);
1337 #elif BOOST_UBLAS_LEGACY_BANDED
1339 size_type l = upper_ + i - j;
1341 l < lower_ + 1 + upper_)
1342 return data () (i, j);
1344 if (is_element_in_band( i, j))
1345 return data () (i, j);
1350 reference operator () (size_type i, size_type j) {
1351 BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1352 BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1353 #ifdef BOOST_UBLAS_OWN_BANDED
1354 size_type k = (std::max) (i, j);
1355 size_type l = lower_ + j - i;
1356 if (k < (std::max) (size1 (), size2 ()) &&
1357 l < lower_ + 1 + upper_)
1358 return data () (i, j);
1359 #elif BOOST_UBLAS_LEGACY_BANDED
1361 size_type l = upper_ + i - j;
1363 l < lower_ + 1 + upper_)
1364 return data () (i, j);
1366 if (is_element_in_band( i, j))
1367 return data () (i, j);
1369 #ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
1370 bad_index ().raise ();
1372 return const_cast<reference>(zero_);
1376 reference operator () (size_type i, size_type j) const {
1377 BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1378 BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1379 #ifdef BOOST_UBLAS_OWN_BANDED
1380 size_type k = (std::max) (i, j);
1381 size_type l = lower_ + j - i;
1382 if (k < (std::max) (size1 (), size2 ()) &&
1383 l < lower_ + 1 + upper_)
1384 return data () (i, j);
1385 #elif BOOST_UBLAS_LEGACY_BANDED
1387 size_type l = upper_ + i - j;
1389 l < lower_ + 1 + upper_)
1390 return data () (i, j);
1392 if (is_element_in_band( i, j))
1393 return data () (i, j);
1395 #ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
1396 bad_index ().raise ();
1398 return const_cast<reference>(zero_);
1404 banded_adaptor &operator = (const banded_adaptor &m) {
1405 matrix_assign<scalar_assign> (*this, m);
1409 banded_adaptor &assign_temporary (banded_adaptor &m) {
1415 banded_adaptor &operator = (const matrix_expression<AE> &ae) {
1416 matrix_assign<scalar_assign> (*this, matrix<value_type> (ae));
1421 banded_adaptor &assign (const matrix_expression<AE> &ae) {
1422 matrix_assign<scalar_assign> (*this, ae);
1427 banded_adaptor& operator += (const matrix_expression<AE> &ae) {
1428 matrix_assign<scalar_assign> (*this, matrix<value_type> (*this + ae));
1433 banded_adaptor &plus_assign (const matrix_expression<AE> &ae) {
1434 matrix_assign<scalar_plus_assign> (*this, ae);
1439 banded_adaptor& operator -= (const matrix_expression<AE> &ae) {
1440 matrix_assign<scalar_assign> (*this, matrix<value_type> (*this - ae));
1445 banded_adaptor &minus_assign (const matrix_expression<AE> &ae) {
1446 matrix_assign<scalar_minus_assign> (*this, ae);
1451 banded_adaptor& operator *= (const AT &at) {
1452 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1457 banded_adaptor& operator /= (const AT &at) {
1458 matrix_assign_scalar<scalar_divides_assign> (*this, at);
1462 // Closure comparison
1464 bool same_closure (const banded_adaptor &ba) const {
1465 return (*this).data ().same_closure (ba.data ());
1470 void swap (banded_adaptor &m) {
1472 BOOST_UBLAS_CHECK (lower_ == m.lower_, bad_size ());
1473 BOOST_UBLAS_CHECK (upper_ == m.upper_, bad_size ());
1474 matrix_swap<scalar_swap> (*this, m);
1478 friend void swap (banded_adaptor &m1, banded_adaptor &m2) {
1484 // Use the matrix iterator
1485 typedef typename M::const_iterator1 const_subiterator1_type;
1486 typedef typename boost::mpl::if_<boost::is_const<M>,
1487 typename M::const_iterator1,
1488 typename M::iterator1>::type subiterator1_type;
1489 typedef typename M::const_iterator2 const_subiterator2_type;
1490 typedef typename boost::mpl::if_<boost::is_const<M>,
1491 typename M::const_iterator2,
1492 typename M::iterator2>::type subiterator2_type;
1495 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1496 typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
1497 typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
1498 typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
1499 typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
1501 class const_iterator1;
1503 class const_iterator2;
1506 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1507 typedef reverse_iterator_base1<iterator1> reverse_iterator1;
1508 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1509 typedef reverse_iterator_base2<iterator2> reverse_iterator2;
1513 const_iterator1 find1 (int rank, size_type i, size_type j) const {
1515 size_type lower_i = (std::max) (difference_type (j - upper_), difference_type (0));
1516 i = (std::max) (i, lower_i);
1517 size_type upper_i = (std::min) (j + 1 + lower_, size1 ());
1518 i = (std::min) (i, upper_i);
1520 return const_iterator1 (*this, data ().find1 (rank, i, j));
1523 iterator1 find1 (int rank, size_type i, size_type j) {
1525 size_type lower_i = (std::max) (difference_type (j - upper_), difference_type (0));
1526 i = (std::max) (i, lower_i);
1527 size_type upper_i = (std::min) (j + 1 + lower_, size1 ());
1528 i = (std::min) (i, upper_i);
1530 return iterator1 (*this, data ().find1 (rank, i, j));
1533 const_iterator2 find2 (int rank, size_type i, size_type j) const {
1535 size_type lower_j = (std::max) (difference_type (i - lower_), difference_type (0));
1536 j = (std::max) (j, lower_j);
1537 size_type upper_j = (std::min) (i + 1 + upper_, size2 ());
1538 j = (std::min) (j, upper_j);
1540 return const_iterator2 (*this, data ().find2 (rank, i, j));
1543 iterator2 find2 (int rank, size_type i, size_type j) {
1545 size_type lower_j = (std::max) (difference_type (i - lower_), difference_type (0));
1546 j = (std::max) (j, lower_j);
1547 size_type upper_j = (std::min) (i + 1 + upper_, size2 ());
1548 j = (std::min) (j, upper_j);
1550 return iterator2 (*this, data ().find2 (rank, i, j));
1553 // Iterators simply are indices.
1555 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1556 class const_iterator1:
1557 public container_const_reference<banded_adaptor>,
1558 public random_access_iterator_base<typename iterator_restrict_traits<
1559 typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1560 const_iterator1, value_type> {
1562 typedef typename const_subiterator1_type::value_type value_type;
1563 typedef typename const_subiterator1_type::difference_type difference_type;
1564 typedef typename const_subiterator1_type::reference reference;
1565 typedef typename const_subiterator1_type::pointer pointer;
1567 typedef const_iterator2 dual_iterator_type;
1568 typedef const_reverse_iterator2 dual_reverse_iterator_type;
1570 // Construction and destruction
1573 container_const_reference<self_type> (), it1_ () {}
1575 const_iterator1 (const self_type &m, const const_subiterator1_type &it1):
1576 container_const_reference<self_type> (m), it1_ (it1) {}
1578 const_iterator1 (const iterator1 &it):
1579 container_const_reference<self_type> (it ()), it1_ (it.it1_) {}
1583 const_iterator1 &operator ++ () {
1588 const_iterator1 &operator -- () {
1593 const_iterator1 &operator += (difference_type n) {
1598 const_iterator1 &operator -= (difference_type n) {
1603 difference_type operator - (const const_iterator1 &it) const {
1604 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1605 return it1_ - it.it1_;
1610 const_reference operator * () const {
1611 size_type i = index1 ();
1612 size_type j = index2 ();
1613 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1614 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1615 #ifdef BOOST_UBLAS_OWN_BANDED
1616 size_type k = (std::max) (i, j);
1617 size_type l = (*this) ().lower () + j - i;
1618 if (k < (std::max) ((*this) ().size1 (), (*this) ().size2 ()) &&
1619 l < (*this) ().lower () + 1 + (*this) ().upper ())
1623 size_type l = (*this) ().upper () + i - j;
1624 if (k < (*this) ().size2 () &&
1625 l < (*this) ().lower () + 1 + (*this) ().upper ())
1628 return (*this) () (i, j);
1631 const_reference operator [] (difference_type n) const {
1632 return *(*this + n);
1635 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1637 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1638 typename self_type::
1640 const_iterator2 begin () const {
1641 return (*this) ().find2 (1, index1 (), 0);
1644 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1645 typename self_type::
1647 const_iterator2 cbegin () const {
1651 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1652 typename self_type::
1654 const_iterator2 end () const {
1655 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1658 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1659 typename self_type::
1661 const_iterator2 cend () const {
1665 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1666 typename self_type::
1668 const_reverse_iterator2 rbegin () const {
1669 return const_reverse_iterator2 (end ());
1672 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1673 typename self_type::
1675 const_reverse_iterator2 crbegin () const {
1679 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1680 typename self_type::
1682 const_reverse_iterator2 rend () const {
1683 return const_reverse_iterator2 (begin ());
1686 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1687 typename self_type::
1689 const_reverse_iterator2 crend () const {
1696 size_type index1 () const {
1697 return it1_.index1 ();
1700 size_type index2 () const {
1701 return it1_.index2 ();
1706 const_iterator1 &operator = (const const_iterator1 &it) {
1707 container_const_reference<self_type>::assign (&it ());
1714 bool operator == (const const_iterator1 &it) const {
1715 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1716 return it1_ == it.it1_;
1719 bool operator < (const const_iterator1 &it) const {
1720 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1721 return it1_ < it.it1_;
1725 const_subiterator1_type it1_;
1730 const_iterator1 begin1 () const {
1731 return find1 (0, 0, 0);
1734 const_iterator1 cbegin1 () const {
1738 const_iterator1 end1 () const {
1739 return find1 (0, size1 (), 0);
1742 const_iterator1 cend1 () const {
1746 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1748 public container_reference<banded_adaptor>,
1749 public random_access_iterator_base<typename iterator_restrict_traits<
1750 typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1751 iterator1, value_type> {
1753 typedef typename subiterator1_type::value_type value_type;
1754 typedef typename subiterator1_type::difference_type difference_type;
1755 typedef typename subiterator1_type::reference reference;
1756 typedef typename subiterator1_type::pointer pointer;
1758 typedef iterator2 dual_iterator_type;
1759 typedef reverse_iterator2 dual_reverse_iterator_type;
1761 // Construction and destruction
1764 container_reference<self_type> (), it1_ () {}
1766 iterator1 (self_type &m, const subiterator1_type &it1):
1767 container_reference<self_type> (m), it1_ (it1) {}
1771 iterator1 &operator ++ () {
1776 iterator1 &operator -- () {
1781 iterator1 &operator += (difference_type n) {
1786 iterator1 &operator -= (difference_type n) {
1791 difference_type operator - (const iterator1 &it) const {
1792 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1793 return it1_ - it.it1_;
1798 reference operator * () const {
1799 size_type i = index1 ();
1800 size_type j = index2 ();
1801 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1802 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1803 #ifdef BOOST_UBLAS_OWN_BANDED
1804 size_type k = (std::max) (i, j);
1805 size_type l = (*this) ().lower () + j - i;
1806 if (k < (std::max) ((*this) ().size1 (), (*this) ().size2 ()) &&
1807 l < (*this) ().lower () + 1 + (*this) ().upper ())
1811 size_type l = (*this) ().upper () + i - j;
1812 if (k < (*this) ().size2 () &&
1813 l < (*this) ().lower () + 1 + (*this) ().upper ())
1816 return (*this) () (i, j);
1819 reference operator [] (difference_type n) const {
1820 return *(*this + n);
1823 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1825 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1826 typename self_type::
1828 iterator2 begin () const {
1829 return (*this) ().find2 (1, index1 (), 0);
1832 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1833 typename self_type::
1835 iterator2 end () const {
1836 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1839 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1840 typename self_type::
1842 reverse_iterator2 rbegin () const {
1843 return reverse_iterator2 (end ());
1846 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1847 typename self_type::
1849 reverse_iterator2 rend () const {
1850 return reverse_iterator2 (begin ());
1856 size_type index1 () const {
1857 return it1_.index1 ();
1860 size_type index2 () const {
1861 return it1_.index2 ();
1866 iterator1 &operator = (const iterator1 &it) {
1867 container_reference<self_type>::assign (&it ());
1874 bool operator == (const iterator1 &it) const {
1875 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1876 return it1_ == it.it1_;
1879 bool operator < (const iterator1 &it) const {
1880 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1881 return it1_ < it.it1_;
1885 subiterator1_type it1_;
1887 friend class const_iterator1;
1892 iterator1 begin1 () {
1893 return find1 (0, 0, 0);
1897 return find1 (0, size1 (), 0);
1900 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1901 class const_iterator2:
1902 public container_const_reference<banded_adaptor>,
1903 public random_access_iterator_base<packed_random_access_iterator_tag,
1904 const_iterator2, value_type> {
1906 typedef typename iterator_restrict_traits<typename const_subiterator2_type::iterator_category,
1907 packed_random_access_iterator_tag>::iterator_category iterator_category;
1908 typedef typename const_subiterator2_type::value_type value_type;
1909 typedef typename const_subiterator2_type::difference_type difference_type;
1910 typedef typename const_subiterator2_type::reference reference;
1911 typedef typename const_subiterator2_type::pointer pointer;
1913 typedef const_iterator1 dual_iterator_type;
1914 typedef const_reverse_iterator1 dual_reverse_iterator_type;
1916 // Construction and destruction
1919 container_const_reference<self_type> (), it2_ () {}
1921 const_iterator2 (const self_type &m, const const_subiterator2_type &it2):
1922 container_const_reference<self_type> (m), it2_ (it2) {}
1924 const_iterator2 (const iterator2 &it):
1925 container_const_reference<self_type> (it ()), it2_ (it.it2_) {}
1929 const_iterator2 &operator ++ () {
1934 const_iterator2 &operator -- () {
1939 const_iterator2 &operator += (difference_type n) {
1944 const_iterator2 &operator -= (difference_type n) {
1949 difference_type operator - (const const_iterator2 &it) const {
1950 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1951 return it2_ - it.it2_;
1956 const_reference operator * () const {
1957 size_type i = index1 ();
1958 size_type j = index2 ();
1959 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1960 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1961 #ifdef BOOST_UBLAS_OWN_BANDED
1962 size_type k = (std::max) (i, j);
1963 size_type l = (*this) ().lower () + j - i;
1964 if (k < (std::max) ((*this) ().size1 (), (*this) ().size2 ()) &&
1965 l < (*this) ().lower () + 1 + (*this) ().upper ())
1969 size_type l = (*this) ().upper () + i - j;
1970 if (k < (*this) ().size2 () &&
1971 l < (*this) ().lower () + 1 + (*this) ().upper ())
1974 return (*this) () (i, j);
1977 const_reference operator [] (difference_type n) const {
1978 return *(*this + n);
1981 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1983 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1984 typename self_type::
1986 const_iterator1 begin () const {
1987 return (*this) ().find1 (1, 0, index2 ());
1990 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1991 typename self_type::
1993 const_iterator1 cbegin () const {
1997 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1998 typename self_type::
2000 const_iterator1 end () const {
2001 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2004 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2005 typename self_type::
2007 const_iterator1 cend () const {
2011 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2012 typename self_type::
2014 const_reverse_iterator1 rbegin () const {
2015 return const_reverse_iterator1 (end ());
2018 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2019 typename self_type::
2021 const_reverse_iterator1 crbegin () const {
2025 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2026 typename self_type::
2028 const_reverse_iterator1 rend () const {
2029 return const_reverse_iterator1 (begin ());
2032 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2033 typename self_type::
2035 const_reverse_iterator1 crend () const {
2042 size_type index1 () const {
2043 return it2_.index1 ();
2046 size_type index2 () const {
2047 return it2_.index2 ();
2052 const_iterator2 &operator = (const const_iterator2 &it) {
2053 container_const_reference<self_type>::assign (&it ());
2060 bool operator == (const const_iterator2 &it) const {
2061 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2062 return it2_ == it.it2_;
2065 bool operator < (const const_iterator2 &it) const {
2066 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2067 return it2_ < it.it2_;
2071 const_subiterator2_type it2_;
2076 const_iterator2 begin2 () const {
2077 return find2 (0, 0, 0);
2080 const_iterator2 cbegin2 () const {
2084 const_iterator2 end2 () const {
2085 return find2 (0, 0, size2 ());
2088 const_iterator2 cend2 () const {
2092 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2094 public container_reference<banded_adaptor>,
2095 public random_access_iterator_base<typename iterator_restrict_traits<
2096 typename subiterator2_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
2097 iterator2, value_type> {
2099 typedef typename subiterator2_type::value_type value_type;
2100 typedef typename subiterator2_type::difference_type difference_type;
2101 typedef typename subiterator2_type::reference reference;
2102 typedef typename subiterator2_type::pointer pointer;
2104 typedef iterator1 dual_iterator_type;
2105 typedef reverse_iterator1 dual_reverse_iterator_type;
2107 // Construction and destruction
2110 container_reference<self_type> (), it2_ () {}
2112 iterator2 (self_type &m, const subiterator2_type &it2):
2113 container_reference<self_type> (m), it2_ (it2) {}
2117 iterator2 &operator ++ () {
2122 iterator2 &operator -- () {
2127 iterator2 &operator += (difference_type n) {
2132 iterator2 &operator -= (difference_type n) {
2137 difference_type operator - (const iterator2 &it) const {
2138 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2139 return it2_ - it.it2_;
2144 reference operator * () const {
2145 size_type i = index1 ();
2146 size_type j = index2 ();
2147 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
2148 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
2149 #ifdef BOOST_UBLAS_OWN_BANDED
2150 size_type k = (std::max) (i, j);
2151 size_type l = (*this) ().lower () + j - i;
2152 if (k < (std::max) ((*this) ().size1 (), (*this) ().size2 ()) &&
2153 l < (*this) ().lower () + 1 + (*this) ().upper ())
2157 size_type l = (*this) ().upper () + i - j;
2158 if (k < (*this) ().size2 () &&
2159 l < (*this) ().lower () + 1 + (*this) ().upper ())
2162 return (*this) () (i, j);
2165 reference operator [] (difference_type n) const {
2166 return *(*this + n);
2169 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2171 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2172 typename self_type::
2174 iterator1 begin () const {
2175 return (*this) ().find1 (1, 0, index2 ());
2178 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2179 typename self_type::
2181 iterator1 end () const {
2182 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2185 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2186 typename self_type::
2188 reverse_iterator1 rbegin () const {
2189 return reverse_iterator1 (end ());
2192 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2193 typename self_type::
2195 reverse_iterator1 rend () const {
2196 return reverse_iterator1 (begin ());
2202 size_type index1 () const {
2203 return it2_.index1 ();
2206 size_type index2 () const {
2207 return it2_.index2 ();
2212 iterator2 &operator = (const iterator2 &it) {
2213 container_reference<self_type>::assign (&it ());
2220 bool operator == (const iterator2 &it) const {
2221 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2222 return it2_ == it.it2_;
2225 bool operator < (const iterator2 &it) const {
2226 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2227 return it2_ < it.it2_;
2231 subiterator2_type it2_;
2233 friend class const_iterator2;
2238 iterator2 begin2 () {
2239 return find2 (0, 0, 0);
2243 return find2 (0, 0, size2 ());
2246 // Reverse iterators
2249 const_reverse_iterator1 rbegin1 () const {
2250 return const_reverse_iterator1 (end1 ());
2253 const_reverse_iterator1 crbegin1 () const {
2257 const_reverse_iterator1 rend1 () const {
2258 return const_reverse_iterator1 (begin1 ());
2261 const_reverse_iterator1 crend1 () const {
2266 reverse_iterator1 rbegin1 () {
2267 return reverse_iterator1 (end1 ());
2270 reverse_iterator1 rend1 () {
2271 return reverse_iterator1 (begin1 ());
2275 const_reverse_iterator2 rbegin2 () const {
2276 return const_reverse_iterator2 (end2 ());
2279 const_reverse_iterator2 crbegin2 () const {
2283 const_reverse_iterator2 rend2 () const {
2284 return const_reverse_iterator2 (begin2 ());
2287 const_reverse_iterator2 crend2 () const {
2292 reverse_iterator2 rbegin2 () {
2293 return reverse_iterator2 (end2 ());
2296 reverse_iterator2 rend2 () {
2297 return reverse_iterator2 (begin2 ());
2301 matrix_closure_type data_;
2304 typedef const value_type const_value_type;
2305 static const_value_type zero_;
2308 // Specialization for temporary_traits
2310 struct vector_temporary_traits< banded_adaptor<M> >
2311 : vector_temporary_traits< M > {} ;
2313 struct vector_temporary_traits< const banded_adaptor<M> >
2314 : vector_temporary_traits< M > {} ;
2317 struct matrix_temporary_traits< banded_adaptor<M> >
2318 : matrix_temporary_traits< M > {} ;
2320 struct matrix_temporary_traits< const banded_adaptor<M> >
2321 : matrix_temporary_traits< M > {} ;
2325 typename banded_adaptor<M>::const_value_type banded_adaptor<M>::zero_ = value_type/*zero*/();
2327 /** \brief A diagonal matrix adaptator: convert a any matrix into a diagonal matrix expression
2329 * For a \f$(m\times m)\f$-dimensional matrix, the \c diagonal_adaptor will provide a diagonal matrix
2330 * with \f$0 \leq i < m\f$ and \f$0 \leq j < m\f$, if \f$i\neq j\f$ then \f$b_{i,j}=0\f$.
2332 * Storage and location are based on those of the underlying matrix. This is important because
2333 * a \c diagonal_adaptor does not copy the matrix data to a new place. Therefore, modifying values
2334 * in a \c diagonal_adaptor matrix will also modify the underlying matrix too.
2336 * \tparam M the type of matrix used to generate the diagonal matrix
2340 class diagonal_adaptor:
2341 public banded_adaptor<M> {
2343 typedef M matrix_type;
2344 typedef banded_adaptor<M> adaptor_type;
2346 // Construction and destruction
2348 diagonal_adaptor ():
2351 diagonal_adaptor (matrix_type &data):
2352 adaptor_type (data) {}
2354 ~diagonal_adaptor () {}
2358 diagonal_adaptor &operator = (const diagonal_adaptor &m) {
2359 adaptor_type::operator = (m);
2364 diagonal_adaptor &operator = (const matrix_expression<AE> &ae) {
2365 adaptor_type::operator = (ae);