Gitiles
Code Review Sign In
realtimeroboticsgroup.org / RealtimeRoboticsGroup / test / dc6866bfe21cff39c0481df3fb8a37d6f44f3e4d / . / include / boost / numeric / ublas / banded.hpp
blob: d135650a983b6cf3556d1e4e697e072f31a2643a [file] [log] [blame]
Brian Silvermandc6866b2018-08-05 00:18:23 -0700[diff] [blame^]1//
2// Copyright (c) 2000-2013
3// Joerg Walter, Mathias Koch, Athanasios Iliopoulos
4//
5// Distributed under the Boost Software License, Version 1.0. (See
6// accompanying file LICENSE_1_0.txt or copy at
7// http://www.boost.org/LICENSE_1_0.txt)
8//
9// The authors gratefully acknowledge the support of
10// GeNeSys mbH & Co. KG in producing this work.
11//
12
13#ifndef _BOOST_UBLAS_BANDED_
14#define _BOOST_UBLAS_BANDED_
15
16#include <boost/numeric/ublas/matrix.hpp>
17#include <boost/numeric/ublas/detail/temporary.hpp>
18
19// Iterators based on ideas of Jeremy Siek
20
21namespace boost { namespace numeric { namespace ublas {
22
23
24namespace hidden {
25
26
27
28/** \brief A helper for band_matrix indexing.
29 *
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;
32 */
33template <class LayoutType>
34class banded_indexing { };
35
36/** \brief A helper for indexing column major banded matrices.
37 *
38 */
39template <>
40class banded_indexing<column_major_tag> {
41public:
42
43 template <class T>
44 BOOST_UBLAS_INLINE static T size(T /*size1*/, T size2) {
45 return size2;
46 }
47
48// template <class T>
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
51// }
52
53 template <class T>
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);
56 }
57};
58
59/** \brief A helper for indexing row major banded matrices.
60 *
61 */
62template <>
63class banded_indexing<row_major_tag> {
64public:
65
66 template <class T>
67 BOOST_UBLAS_INLINE static T size(T size1, T /*size2*/) {
68 return size1;
69 }
70
71 // template <class T>
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
74 // }
75
76 template <class T>
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);
79 }
80};
81
82}
83
84 /** \brief A banded matrix of values of type \c T.
85 *
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.
91 *
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
95 */
96 template<class T, class L, class A>
97 class banded_matrix:
98 public matrix_container<banded_matrix<T, L, A> > {
99
100 typedef T *pointer;
101 typedef L layout_type;
102 typedef banded_matrix<T, L, A> self_type;
103
104
105
106 public:
107#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
108 using matrix_container<self_type>::operator ();
109#endif
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;
122
123 private:
124 public:
125
126 // Construction and destruction
127 BOOST_UBLAS_INLINE
128 banded_matrix ():
129 matrix_container<self_type> (),
130 size1_ (0), size2_ (0),
131 lower_ (0), upper_ (0), data_ (0) {}
132 BOOST_UBLAS_INLINE
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))
139#else
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
141#endif
142 {
143 }
144 BOOST_UBLAS_INLINE
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) {}
149 BOOST_UBLAS_INLINE
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_) {}
154 template<class AE>
155 BOOST_UBLAS_INLINE
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_))
162#else
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
164#endif
165 {
166 matrix_assign<scalar_assign> (*this, ae);
167 }
168
169 // Accessors
170 BOOST_UBLAS_INLINE
171 size_type size1 () const {
172 return size1_;
173 }
174 BOOST_UBLAS_INLINE
175 size_type size2 () const {
176 return size2_;
177 }
178 BOOST_UBLAS_INLINE
179 size_type lower () const {
180 return lower_;
181 }
182 BOOST_UBLAS_INLINE
183 size_type upper () const {
184 return upper_;
185 }
186
187 // Storage accessors
188 BOOST_UBLAS_INLINE
189 const array_type &data () const {
190 return data_;
191 }
192 BOOST_UBLAS_INLINE
193 array_type &data () {
194 return data_;
195 }
196
197#if !defined (BOOST_UBLAS_OWN_BANDED)||(BOOST_UBLAS_LEGACY_BANDED)
198 BOOST_UBLAS_INLINE
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
202 }
203#endif
204 // Resizing
205 BOOST_UBLAS_INLINE
206 void resize (size_type size1, size_type size2, size_type lower = 0, size_type upper = 0, bool preserve = true) {
207 if (preserve) {
208 self_type temporary (size1, size2, lower, upper);
209 detail::matrix_resize_preserve<layout_type> (*this, temporary);
210 }
211 else {
212 data ().resize ((std::max) (size1, size2) * (lower + 1 + upper));
213 size1_ = size1;
214 size2_ = size2;
215 lower_ = lower;
216 upper_ = upper;
217 }
218 }
219
220 BOOST_UBLAS_INLINE
221 void resize_packed_preserve (size_type size1, size_type size2, size_type lower = 0, size_type upper = 0) {
222 size1_ = size1;
223 size2_ = size2;
224 lower_ = lower;
225 upper_ = upper;
226 data ().resize ((std::max) (size1, size2) * (lower + 1 + upper), value_type ());
227 }
228
229 // Element access
230 BOOST_UBLAS_INLINE
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;
244 if (k < size2_ &&
245 l < lower_ + 1 + upper_)
246 return data () [layout_type::element (k, size2_,
247 l, lower_ + 1 + upper_)];
248#else // New default
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)];
252 }
253#endif
254 return zero_;
255 }
256
257 BOOST_UBLAS_INLINE
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;
269 if (! (k < size2_ &&
270 l < lower_ + 1 + upper_) ) {
271 bad_index ().raise ();
272 // NEVER reached
273 }
274 return data () [layout_type::element (k, size2_,
275 l, lower_ + 1 + upper_)];
276#else
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)];
280#endif
281 }
282 BOOST_UBLAS_INLINE
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 ();
292 // NEVER reached
293 }
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;
299 if (! (k < size2_ &&
300 l < lower_ + 1 + upper_) ) {
301 bad_index ().raise ();
302 // NEVER reached
303 }
304 return data () [layout_type::element (k, size2_,
305 l, lower_ + 1 + upper_)];
306#else
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)];
310#endif
311
312 }
313
314 // Element assignment
315 BOOST_UBLAS_INLINE
316 reference insert_element (size_type i, size_type j, const_reference t) {
317 return (operator () (i, j) = t);
318 }
319 BOOST_UBLAS_INLINE
320 void erase_element (size_type i, size_type j) {
321 operator () (i, j) = value_type/*zero*/();
322 }
323
324 // Zeroing
325 BOOST_UBLAS_INLINE
326 void clear () {
327 std::fill (data ().begin (), data ().end (), value_type/*zero*/());
328 }
329
330 // Assignment
331 BOOST_UBLAS_INLINE
332 banded_matrix &operator = (const banded_matrix &m) {
333 size1_ = m.size1_;
334 size2_ = m.size2_;
335 lower_ = m.lower_;
336 upper_ = m.upper_;
337 data () = m.data ();
338 return *this;
339 }
340 BOOST_UBLAS_INLINE
341 banded_matrix &assign_temporary (banded_matrix &m) {
342 swap (m);
343 return *this;
344 }
345 template<class AE>
346 BOOST_UBLAS_INLINE
347 banded_matrix &operator = (const matrix_expression<AE> &ae) {
348 self_type temporary (ae, lower_, upper_);
349 return assign_temporary (temporary);
350 }
351 template<class AE>
352 BOOST_UBLAS_INLINE
353 banded_matrix &assign (const matrix_expression<AE> &ae) {
354 matrix_assign<scalar_assign> (*this, ae);
355 return *this;
356 }
357 template<class AE>
358 BOOST_UBLAS_INLINE
359 banded_matrix& operator += (const matrix_expression<AE> &ae) {
360 self_type temporary (*this + ae, lower_, upper_);
361 return assign_temporary (temporary);
362 }
363 template<class AE>
364 BOOST_UBLAS_INLINE
365 banded_matrix &plus_assign (const matrix_expression<AE> &ae) {
366 matrix_assign<scalar_plus_assign> (*this, ae);
367 return *this;
368 }
369 template<class AE>
370 BOOST_UBLAS_INLINE
371 banded_matrix& operator -= (const matrix_expression<AE> &ae) {
372 self_type temporary (*this - ae, lower_, upper_);
373 return assign_temporary (temporary);
374 }
375 template<class AE>
376 BOOST_UBLAS_INLINE
377 banded_matrix &minus_assign (const matrix_expression<AE> &ae) {
378 matrix_assign<scalar_minus_assign> (*this, ae);
379 return *this;
380 }
381 template<class AT>
382 BOOST_UBLAS_INLINE
383 banded_matrix& operator *= (const AT &at) {
384 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
385 return *this;
386 }
387 template<class AT>
388 BOOST_UBLAS_INLINE
389 banded_matrix& operator /= (const AT &at) {
390 matrix_assign_scalar<scalar_divides_assign> (*this, at);
391 return *this;
392 }
393
394 // Swapping
395 BOOST_UBLAS_INLINE
396 void swap (banded_matrix &m) {
397 if (this != &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 ());
403 }
404 }
405 BOOST_UBLAS_INLINE
406 friend void swap (banded_matrix &m1, banded_matrix &m2) {
407 m1.swap (m2);
408 }
409
410 // Iterator types
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;
416#else
417 class const_iterator1;
418 class iterator1;
419 class const_iterator2;
420 class iterator2;
421#endif
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;
426
427 // Element lookup
428 BOOST_UBLAS_INLINE
429 const_iterator1 find1 (int rank, size_type i, size_type j) const {
430 if (rank == 1) {
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);
435 }
436 return const_iterator1 (*this, i, j);
437 }
438 BOOST_UBLAS_INLINE
439 iterator1 find1 (int rank, size_type i, size_type j) {
440 if (rank == 1) {
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);
445 }
446 return iterator1 (*this, i, j);
447 }
448 BOOST_UBLAS_INLINE
449 const_iterator2 find2 (int rank, size_type i, size_type j) const {
450 if (rank == 1) {
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);
455 }
456 return const_iterator2 (*this, i, j);
457 }
458 BOOST_UBLAS_INLINE
459 iterator2 find2 (int rank, size_type i, size_type j) {
460 if (rank == 1) {
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);
465 }
466 return iterator2 (*this, i, j);
467 }
468
469 // Iterators simply are indices.
470
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> {
476 public:
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;
481
482 typedef const_iterator2 dual_iterator_type;
483 typedef const_reverse_iterator2 dual_reverse_iterator_type;
484
485 // Construction and destruction
486 BOOST_UBLAS_INLINE
487 const_iterator1 ():
488 container_const_reference<self_type> (), it1_ (), it2_ () {}
489 BOOST_UBLAS_INLINE
490 const_iterator1 (const self_type &m, size_type it1, size_type it2):
491 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
492 BOOST_UBLAS_INLINE
493 const_iterator1 (const iterator1 &it):
494 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
495
496 // Arithmetic
497 BOOST_UBLAS_INLINE
498 const_iterator1 &operator ++ () {
499 ++ it1_;
500 return *this;
501 }
502 BOOST_UBLAS_INLINE
503 const_iterator1 &operator -- () {
504 -- it1_;
505 return *this;
506 }
507 BOOST_UBLAS_INLINE
508 const_iterator1 &operator += (difference_type n) {
509 it1_ += n;
510 return *this;
511 }
512 BOOST_UBLAS_INLINE
513 const_iterator1 &operator -= (difference_type n) {
514 it1_ -= n;
515 return *this;
516 }
517 BOOST_UBLAS_INLINE
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_;
522 }
523
524 // Dereference
525 BOOST_UBLAS_INLINE
526 const_reference operator * () const {
527 return (*this) () (it1_, it2_);
528 }
529 BOOST_UBLAS_INLINE
530 const_reference operator [] (difference_type n) const {
531 return *(*this + n);
532 }
533
534#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
535 BOOST_UBLAS_INLINE
536#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
537 typename self_type::
538#endif
539 const_iterator2 begin () const {
540 return (*this) ().find2 (1, it1_, 0);
541 }
542 BOOST_UBLAS_INLINE
543#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
544 typename self_type::
545#endif
546 const_iterator2 cbegin () const {
547 return begin ();
548 }
549 BOOST_UBLAS_INLINE
550#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
551 typename self_type::
552#endif
553 const_iterator2 end () const {
554 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
555 }
556 BOOST_UBLAS_INLINE
557#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
558 typename self_type::
559#endif
560 const_iterator2 cend () const {
561 return end ();
562 }
563 BOOST_UBLAS_INLINE
564#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
565 typename self_type::
566#endif
567 const_reverse_iterator2 rbegin () const {
568 return const_reverse_iterator2 (end ());
569 }
570 BOOST_UBLAS_INLINE
571#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
572 typename self_type::
573#endif
574 const_reverse_iterator2 crbegin () const {
575 return rbegin ();
576 }
577 BOOST_UBLAS_INLINE
578#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
579 typename self_type::
580#endif
581 const_reverse_iterator2 rend () const {
582 return const_reverse_iterator2 (begin ());
583 }
584 BOOST_UBLAS_INLINE
585#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
586 typename self_type::
587#endif
588 const_reverse_iterator2 crend () const {
589 return rend ();
590 }
591#endif
592
593 // Indices
594 BOOST_UBLAS_INLINE
595 size_type index1 () const {
596 return it1_;
597 }
598 BOOST_UBLAS_INLINE
599 size_type index2 () const {
600 return it2_;
601 }
602
603 // Assignment
604 BOOST_UBLAS_INLINE
605 const_iterator1 &operator = (const const_iterator1 &it) {
606 container_const_reference<self_type>::assign (&it ());
607 it1_ = it.it1_;
608 it2_ = it.it2_;
609 return *this;
610 }
611
612 // Comparison
613 BOOST_UBLAS_INLINE
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_;
618 }
619 BOOST_UBLAS_INLINE
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_;
624 }
625
626 private:
627 size_type it1_;
628 size_type it2_;
629 };
630#endif
631
632 BOOST_UBLAS_INLINE
633 const_iterator1 begin1 () const {
634 return find1 (0, 0, 0);
635 }
636 BOOST_UBLAS_INLINE
637 const_iterator1 cbegin1 () const {
638 return begin1 ();
639 }
640 BOOST_UBLAS_INLINE
641 const_iterator1 end1 () const {
642 return find1 (0, size1_, 0);
643 }
644 BOOST_UBLAS_INLINE
645 const_iterator1 cend1 () const {
646 return end1 ();
647 }
648
649#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
650 class iterator1:
651 public container_reference<banded_matrix>,
652 public random_access_iterator_base<packed_random_access_iterator_tag,
653 iterator1, value_type> {
654 public:
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;
659
660 typedef iterator2 dual_iterator_type;
661 typedef reverse_iterator2 dual_reverse_iterator_type;
662
663 // Construction and destruction
664 BOOST_UBLAS_INLINE
665 iterator1 ():
666 container_reference<self_type> (), it1_ (), it2_ () {}
667 BOOST_UBLAS_INLINE
668 iterator1 (self_type &m, size_type it1, size_type it2):
669 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
670
671 // Arithmetic
672 BOOST_UBLAS_INLINE
673 iterator1 &operator ++ () {
674 ++ it1_;
675 return *this;
676 }
677 BOOST_UBLAS_INLINE
678 iterator1 &operator -- () {
679 -- it1_;
680 return *this;
681 }
682 BOOST_UBLAS_INLINE
683 iterator1 &operator += (difference_type n) {
684 it1_ += n;
685 return *this;
686 }
687 BOOST_UBLAS_INLINE
688 iterator1 &operator -= (difference_type n) {
689 it1_ -= n;
690 return *this;
691 }
692 BOOST_UBLAS_INLINE
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_;
697 }
698
699 // Dereference
700 BOOST_UBLAS_INLINE
701 reference operator * () const {
702 return (*this) ().at_element (it1_, it2_);
703 }
704 BOOST_UBLAS_INLINE
705 reference operator [] (difference_type n) const {
706 return *(*this + n);
707 }
708
709#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
710 BOOST_UBLAS_INLINE
711#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
712 typename self_type::
713#endif
714 iterator2 begin () const {
715 return (*this) ().find2 (1, it1_, 0);
716 }
717 BOOST_UBLAS_INLINE
718#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
719 typename self_type::
720#endif
721 iterator2 end () const {
722 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
723 }
724
725 BOOST_UBLAS_INLINE
726#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
727 typename self_type::
728#endif
729 reverse_iterator2 rbegin () const {
730 return reverse_iterator2 (end ());
731 }
732 BOOST_UBLAS_INLINE
733#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
734 typename self_type::
735#endif
736 reverse_iterator2 rend () const {
737 return reverse_iterator2 (begin ());
738 }
739#endif
740
741 // Indices
742 BOOST_UBLAS_INLINE
743 size_type index1 () const {
744 return it1_;
745 }
746 BOOST_UBLAS_INLINE
747 size_type index2 () const {
748 return it2_;
749 }
750
751 // Assignment
752 BOOST_UBLAS_INLINE
753 iterator1 &operator = (const iterator1 &it) {
754 container_reference<self_type>::assign (&it ());
755 it1_ = it.it1_;
756 it2_ = it.it2_;
757 return *this;
758 }
759
760 // Comparison
761 BOOST_UBLAS_INLINE
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_;
766 }
767 BOOST_UBLAS_INLINE
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_;
772 }
773
774 private:
775 size_type it1_;
776 size_type it2_;
777
778 friend class const_iterator1;
779 };
780#endif
781
782 BOOST_UBLAS_INLINE
783 iterator1 begin1 () {
784 return find1 (0, 0, 0);
785 }
786 BOOST_UBLAS_INLINE
787 iterator1 end1 () {
788 return find1 (0, size1_, 0);
789 }
790
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> {
796 public:
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;
801
802 typedef const_iterator1 dual_iterator_type;
803 typedef const_reverse_iterator1 dual_reverse_iterator_type;
804
805 // Construction and destruction
806 BOOST_UBLAS_INLINE
807 const_iterator2 ():
808 container_const_reference<self_type> (), it1_ (), it2_ () {}
809 BOOST_UBLAS_INLINE
810 const_iterator2 (const self_type &m, size_type it1, size_type it2):
811 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
812 BOOST_UBLAS_INLINE
813 const_iterator2 (const iterator2 &it):
814 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
815
816 // Arithmetic
817 BOOST_UBLAS_INLINE
818 const_iterator2 &operator ++ () {
819 ++ it2_;
820 return *this;
821 }
822 BOOST_UBLAS_INLINE
823 const_iterator2 &operator -- () {
824 -- it2_;
825 return *this;
826 }
827 BOOST_UBLAS_INLINE
828 const_iterator2 &operator += (difference_type n) {
829 it2_ += n;
830 return *this;
831 }
832 BOOST_UBLAS_INLINE
833 const_iterator2 &operator -= (difference_type n) {
834 it2_ -= n;
835 return *this;
836 }
837 BOOST_UBLAS_INLINE
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_;
842 }
843
844 // Dereference
845 BOOST_UBLAS_INLINE
846 const_reference operator * () const {
847 return (*this) () (it1_, it2_);
848 }
849 BOOST_UBLAS_INLINE
850 const_reference operator [] (difference_type n) const {
851 return *(*this + n);
852 }
853
854#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
855 BOOST_UBLAS_INLINE
856#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
857 typename self_type::
858#endif
859 const_iterator1 begin () const {
860 return (*this) ().find1 (1, 0, it2_);
861 }
862 BOOST_UBLAS_INLINE
863#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
864 typename self_type::
865#endif
866 const_iterator1 cbegin () const {
867 return begin ();
868 }
869 BOOST_UBLAS_INLINE
870#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
871 typename self_type::
872#endif
873 const_iterator1 end () const {
874 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
875 }
876 BOOST_UBLAS_INLINE
877#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
878 typename self_type::
879#endif
880 const_iterator1 cend () const {
881 return end();
882 }
883
884 BOOST_UBLAS_INLINE
885#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
886 typename self_type::
887#endif
888 const_reverse_iterator1 rbegin () const {
889 return const_reverse_iterator1 (end ());
890 }
891 BOOST_UBLAS_INLINE
892#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
893 typename self_type::
894#endif
895 const_reverse_iterator1 crbegin () const {
896 return rbegin ();
897 }
898 BOOST_UBLAS_INLINE
899#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
900 typename self_type::
901#endif
902 const_reverse_iterator1 rend () const {
903 return const_reverse_iterator1 (begin ());
904 }
905 BOOST_UBLAS_INLINE
906#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
907 typename self_type::
908#endif
909 const_reverse_iterator1 crend () const {
910 return rend ();
911 }
912#endif
913
914 // Indices
915 BOOST_UBLAS_INLINE
916 size_type index1 () const {
917 return it1_;
918 }
919 BOOST_UBLAS_INLINE
920 size_type index2 () const {
921 return it2_;
922 }
923
924 // Assignment
925 BOOST_UBLAS_INLINE
926 const_iterator2 &operator = (const const_iterator2 &it) {
927 container_const_reference<self_type>::assign (&it ());
928 it1_ = it.it1_;
929 it2_ = it.it2_;
930 return *this;
931 }
932
933 // Comparison
934 BOOST_UBLAS_INLINE
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_;
939 }
940 BOOST_UBLAS_INLINE
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_;
945 }
946
947 private:
948 size_type it1_;
949 size_type it2_;
950 };
951#endif
952
953 BOOST_UBLAS_INLINE
954 const_iterator2 begin2 () const {
955 return find2 (0, 0, 0);
956 }
957 BOOST_UBLAS_INLINE
958 const_iterator2 cbegin2 () const {
959 return begin2 ();
960 }
961 BOOST_UBLAS_INLINE
962 const_iterator2 end2 () const {
963 return find2 (0, 0, size2_);
964 }
965 BOOST_UBLAS_INLINE
966 const_iterator2 cend2 () const {
967 return end2 ();
968 }
969
970#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
971 class iterator2:
972 public container_reference<banded_matrix>,
973 public random_access_iterator_base<packed_random_access_iterator_tag,
974 iterator2, value_type> {
975 public:
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;
980
981 typedef iterator1 dual_iterator_type;
982 typedef reverse_iterator1 dual_reverse_iterator_type;
983
984 // Construction and destruction
985 BOOST_UBLAS_INLINE
986 iterator2 ():
987 container_reference<self_type> (), it1_ (), it2_ () {}
988 BOOST_UBLAS_INLINE
989 iterator2 (self_type &m, size_type it1, size_type it2):
990 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
991
992 // Arithmetic
993 BOOST_UBLAS_INLINE
994 iterator2 &operator ++ () {
995 ++ it2_;
996 return *this;
997 }
998 BOOST_UBLAS_INLINE
999 iterator2 &operator -- () {
1000 -- it2_;
1001 return *this;
1002 }
1003 BOOST_UBLAS_INLINE
1004 iterator2 &operator += (difference_type n) {
1005 it2_ += n;
1006 return *this;
1007 }
1008 BOOST_UBLAS_INLINE
1009 iterator2 &operator -= (difference_type n) {
1010 it2_ -= n;
1011 return *this;
1012 }
1013 BOOST_UBLAS_INLINE
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_;
1018 }
1019
1020 // Dereference
1021 BOOST_UBLAS_INLINE
1022 reference operator * () const {
1023 return (*this) ().at_element (it1_, it2_);
1024 }
1025 BOOST_UBLAS_INLINE
1026 reference operator [] (difference_type n) const {
1027 return *(*this + n);
1028 }
1029
1030#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1031 BOOST_UBLAS_INLINE
1032#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1033 typename self_type::
1034#endif
1035 iterator1 begin () const {
1036 return (*this) ().find1 (1, 0, it2_);
1037 }
1038 BOOST_UBLAS_INLINE
1039#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1040 typename self_type::
1041#endif
1042 iterator1 end () const {
1043 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
1044 }
1045 BOOST_UBLAS_INLINE
1046#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1047 typename self_type::
1048#endif
1049 reverse_iterator1 rbegin () const {
1050 return reverse_iterator1 (end ());
1051 }
1052 BOOST_UBLAS_INLINE
1053#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1054 typename self_type::
1055#endif
1056 reverse_iterator1 rend () const {
1057 return reverse_iterator1 (begin ());
1058 }
1059#endif
1060
1061 // Indices
1062 BOOST_UBLAS_INLINE
1063 size_type index1 () const {
1064 return it1_;
1065 }
1066 BOOST_UBLAS_INLINE
1067 size_type index2 () const {
1068 return it2_;
1069 }
1070
1071 // Assignment
1072 BOOST_UBLAS_INLINE
1073 iterator2 &operator = (const iterator2 &it) {
1074 container_reference<self_type>::assign (&it ());
1075 it1_ = it.it1_;
1076 it2_ = it.it2_;
1077 return *this;
1078 }
1079
1080 // Comparison
1081 BOOST_UBLAS_INLINE
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_;
1086 }
1087 BOOST_UBLAS_INLINE
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_;
1092 }
1093
1094 private:
1095 size_type it1_;
1096 size_type it2_;
1097
1098 friend class const_iterator2;
1099 };
1100#endif
1101
1102 BOOST_UBLAS_INLINE
1103 iterator2 begin2 () {
1104 return find2 (0, 0, 0);
1105 }
1106 BOOST_UBLAS_INLINE
1107 iterator2 end2 () {
1108 return find2 (0, 0, size2_);
1109 }
1110
1111 // Reverse iterators
1112
1113 BOOST_UBLAS_INLINE
1114 const_reverse_iterator1 rbegin1 () const {
1115 return const_reverse_iterator1 (end1 ());
1116 }
1117 BOOST_UBLAS_INLINE
1118 const_reverse_iterator1 crbegin1 () const {
1119 return rbegin1 ();
1120 }
1121 BOOST_UBLAS_INLINE
1122 const_reverse_iterator1 rend1 () const {
1123 return const_reverse_iterator1 (begin1 ());
1124 }
1125 BOOST_UBLAS_INLINE
1126 const_reverse_iterator1 crend1 () const {
1127 return rend1 ();
1128 }
1129
1130 BOOST_UBLAS_INLINE
1131 reverse_iterator1 rbegin1 () {
1132 return reverse_iterator1 (end1 ());
1133 }
1134 BOOST_UBLAS_INLINE
1135 reverse_iterator1 rend1 () {
1136 return reverse_iterator1 (begin1 ());
1137 }
1138
1139 BOOST_UBLAS_INLINE
1140 const_reverse_iterator2 rbegin2 () const {
1141 return const_reverse_iterator2 (end2 ());
1142 }
1143 BOOST_UBLAS_INLINE
1144 const_reverse_iterator2 crbegin2 () const {
1145 return rbegin2 ();
1146 }
1147 BOOST_UBLAS_INLINE
1148 const_reverse_iterator2 rend2 () const {
1149 return const_reverse_iterator2 (begin2 ());
1150 }
1151 BOOST_UBLAS_INLINE
1152 const_reverse_iterator2 crend2 () const {
1153 return rend2 ();
1154 }
1155
1156 BOOST_UBLAS_INLINE
1157 reverse_iterator2 rbegin2 () {
1158 return reverse_iterator2 (end2 ());
1159 }
1160 BOOST_UBLAS_INLINE
1161 reverse_iterator2 rend2 () {
1162 return reverse_iterator2 (begin2 ());
1163 }
1164
1165 private:
1166 size_type size1_;
1167 size_type size2_;
1168 size_type lower_;
1169 size_type upper_;
1170 array_type data_;
1171 typedef const value_type const_value_type;
1172 static const_value_type zero_;
1173 };
1174
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*/();
1177
1178
1179 /** \brief A diagonal matrix of values of type \c T, which is a specialization of a banded matrix
1180 *
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.
1184 *
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.
1189 *
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
1193 */
1194 template<class T, class L, class A>
1195 class diagonal_matrix:
1196 public banded_matrix<T, L, A> {
1197 public:
1198 typedef typename A::size_type size_type;
1199 typedef banded_matrix<T, L, A> matrix_type;
1200 typedef A array_type;
1201
1202 // Construction and destruction
1203 BOOST_UBLAS_INLINE
1204 diagonal_matrix ():
1205 matrix_type () {}
1206 BOOST_UBLAS_INLINE
1207 diagonal_matrix (size_type size):
1208 matrix_type (size, size) {}
1209 BOOST_UBLAS_INLINE
1210 diagonal_matrix (size_type size, const array_type& data):
1211 matrix_type (size, size, 0, 0, data) {}
1212 BOOST_UBLAS_INLINE
1213 diagonal_matrix (size_type size1, size_type size2):
1214 matrix_type (size1, size2) {}
1215 template<class AE>
1216 BOOST_UBLAS_INLINE
1217 diagonal_matrix (const matrix_expression<AE> &ae):
1218 matrix_type (ae) {}
1219 BOOST_UBLAS_INLINE
1220 ~diagonal_matrix () {}
1221
1222 // Assignment
1223 BOOST_UBLAS_INLINE
1224 diagonal_matrix &operator = (const diagonal_matrix &m) {
1225 matrix_type::operator = (m);
1226 return *this;
1227 }
1228 template<class AE>
1229 BOOST_UBLAS_INLINE
1230 diagonal_matrix &operator = (const matrix_expression<AE> &ae) {
1231 matrix_type::operator = (ae);
1232 return *this;
1233 }
1234 };
1235
1236 /** \brief A banded matrix adaptator: convert a any matrix into a banded matrix expression
1237 *
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$.
1241 *
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.
1245 *
1246 * \tparam M the type of matrix used to generate a banded matrix
1247 */
1248 template<class M>
1249 class banded_adaptor:
1250 public matrix_expression<banded_adaptor<M> > {
1251
1252 typedef banded_adaptor<M> self_type;
1253
1254 public:
1255#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1256 using matrix_expression<self_type>::operator ();
1257#endif
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;
1278
1279 // Construction and destruction
1280 BOOST_UBLAS_INLINE
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) {}
1284 BOOST_UBLAS_INLINE
1285 banded_adaptor (const banded_adaptor &m):
1286 matrix_expression<self_type> (),
1287 data_ (m.data_), lower_ (m.lower_), upper_ (m.upper_) {}
1288
1289 // Accessors
1290 BOOST_UBLAS_INLINE
1291 size_type size1 () const {
1292 return data_.size1 ();
1293 }
1294 BOOST_UBLAS_INLINE
1295 size_type size2 () const {
1296 return data_.size2 ();
1297 }
1298 BOOST_UBLAS_INLINE
1299 size_type lower () const {
1300 return lower_;
1301 }
1302 BOOST_UBLAS_INLINE
1303 size_type upper () const {
1304 return upper_;
1305 }
1306
1307 // Storage accessors
1308 BOOST_UBLAS_INLINE
1309 const matrix_closure_type &data () const {
1310 return data_;
1311 }
1312 BOOST_UBLAS_INLINE
1313 matrix_closure_type &data () {
1314 return data_;
1315 }
1316
1317#if !defined (BOOST_UBLAS_OWN_BANDED)||(BOOST_UBLAS_LEGACY_BANDED)
1318 BOOST_UBLAS_INLINE
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
1322 }
1323#endif
1324
1325 // Element access
1326#ifndef BOOST_UBLAS_PROXY_CONST_MEMBER
1327 BOOST_UBLAS_INLINE
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
1338 size_type k = j;
1339 size_type l = upper_ + i - j;
1340 if (k < size2 () &&
1341 l < lower_ + 1 + upper_)
1342 return data () (i, j);
1343#else
1344 if (is_element_in_band( i, j))
1345 return data () (i, j);
1346#endif
1347 return zero_;
1348 }
1349 BOOST_UBLAS_INLINE
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
1360 size_type k = j;
1361 size_type l = upper_ + i - j;
1362 if (k < size2 () &&
1363 l < lower_ + 1 + upper_)
1364 return data () (i, j);
1365#else
1366 if (is_element_in_band( i, j))
1367 return data () (i, j);
1368#endif
1369#ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
1370 bad_index ().raise ();
1371#endif
1372 return const_cast<reference>(zero_);
1373 }
1374#else
1375 BOOST_UBLAS_INLINE
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
1386 size_type k = j;
1387 size_type l = upper_ + i - j;
1388 if (k < size2 () &&
1389 l < lower_ + 1 + upper_)
1390 return data () (i, j);
1391#else
1392 if (is_element_in_band( i, j))
1393 return data () (i, j);
1394#endif
1395#ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
1396 bad_index ().raise ();
1397#endif
1398 return const_cast<reference>(zero_);
1399 }
1400#endif
1401
1402 // Assignment
1403 BOOST_UBLAS_INLINE
1404 banded_adaptor &operator = (const banded_adaptor &m) {
1405 matrix_assign<scalar_assign> (*this, m);
1406 return *this;
1407 }
1408 BOOST_UBLAS_INLINE
1409 banded_adaptor &assign_temporary (banded_adaptor &m) {
1410 *this = m;
1411 return *this;
1412 }
1413 template<class AE>
1414 BOOST_UBLAS_INLINE
1415 banded_adaptor &operator = (const matrix_expression<AE> &ae) {
1416 matrix_assign<scalar_assign> (*this, matrix<value_type> (ae));
1417 return *this;
1418 }
1419 template<class AE>
1420 BOOST_UBLAS_INLINE
1421 banded_adaptor &assign (const matrix_expression<AE> &ae) {
1422 matrix_assign<scalar_assign> (*this, ae);
1423 return *this;
1424 }
1425 template<class AE>
1426 BOOST_UBLAS_INLINE
1427 banded_adaptor& operator += (const matrix_expression<AE> &ae) {
1428 matrix_assign<scalar_assign> (*this, matrix<value_type> (*this + ae));
1429 return *this;
1430 }
1431 template<class AE>
1432 BOOST_UBLAS_INLINE
1433 banded_adaptor &plus_assign (const matrix_expression<AE> &ae) {
1434 matrix_assign<scalar_plus_assign> (*this, ae);
1435 return *this;
1436 }
1437 template<class AE>
1438 BOOST_UBLAS_INLINE
1439 banded_adaptor& operator -= (const matrix_expression<AE> &ae) {
1440 matrix_assign<scalar_assign> (*this, matrix<value_type> (*this - ae));
1441 return *this;
1442 }
1443 template<class AE>
1444 BOOST_UBLAS_INLINE
1445 banded_adaptor &minus_assign (const matrix_expression<AE> &ae) {
1446 matrix_assign<scalar_minus_assign> (*this, ae);
1447 return *this;
1448 }
1449 template<class AT>
1450 BOOST_UBLAS_INLINE
1451 banded_adaptor& operator *= (const AT &at) {
1452 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1453 return *this;
1454 }
1455 template<class AT>
1456 BOOST_UBLAS_INLINE
1457 banded_adaptor& operator /= (const AT &at) {
1458 matrix_assign_scalar<scalar_divides_assign> (*this, at);
1459 return *this;
1460 }
1461
1462 // Closure comparison
1463 BOOST_UBLAS_INLINE
1464 bool same_closure (const banded_adaptor &ba) const {
1465 return (*this).data ().same_closure (ba.data ());
1466 }
1467
1468 // Swapping
1469 BOOST_UBLAS_INLINE
1470 void swap (banded_adaptor &m) {
1471 if (this != &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);
1475 }
1476 }
1477 BOOST_UBLAS_INLINE
1478 friend void swap (banded_adaptor &m1, banded_adaptor &m2) {
1479 m1.swap (m2);
1480 }
1481
1482 // Iterator types
1483 private:
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;
1493
1494 public:
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;
1500#else
1501 class const_iterator1;
1502 class iterator1;
1503 class const_iterator2;
1504 class iterator2;
1505#endif
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;
1510
1511 // Element lookup
1512 BOOST_UBLAS_INLINE
1513 const_iterator1 find1 (int rank, size_type i, size_type j) const {
1514 if (rank == 1) {
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);
1519 }
1520 return const_iterator1 (*this, data ().find1 (rank, i, j));
1521 }
1522 BOOST_UBLAS_INLINE
1523 iterator1 find1 (int rank, size_type i, size_type j) {
1524 if (rank == 1) {
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);
1529 }
1530 return iterator1 (*this, data ().find1 (rank, i, j));
1531 }
1532 BOOST_UBLAS_INLINE
1533 const_iterator2 find2 (int rank, size_type i, size_type j) const {
1534 if (rank == 1) {
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);
1539 }
1540 return const_iterator2 (*this, data ().find2 (rank, i, j));
1541 }
1542 BOOST_UBLAS_INLINE
1543 iterator2 find2 (int rank, size_type i, size_type j) {
1544 if (rank == 1) {
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);
1549 }
1550 return iterator2 (*this, data ().find2 (rank, i, j));
1551 }
1552
1553 // Iterators simply are indices.
1554
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> {
1561 public:
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;
1566
1567 typedef const_iterator2 dual_iterator_type;
1568 typedef const_reverse_iterator2 dual_reverse_iterator_type;
1569
1570 // Construction and destruction
1571 BOOST_UBLAS_INLINE
1572 const_iterator1 ():
1573 container_const_reference<self_type> (), it1_ () {}
1574 BOOST_UBLAS_INLINE
1575 const_iterator1 (const self_type &m, const const_subiterator1_type &it1):
1576 container_const_reference<self_type> (m), it1_ (it1) {}
1577 BOOST_UBLAS_INLINE
1578 const_iterator1 (const iterator1 &it):
1579 container_const_reference<self_type> (it ()), it1_ (it.it1_) {}
1580
1581 // Arithmetic
1582 BOOST_UBLAS_INLINE
1583 const_iterator1 &operator ++ () {
1584 ++ it1_;
1585 return *this;
1586 }
1587 BOOST_UBLAS_INLINE
1588 const_iterator1 &operator -- () {
1589 -- it1_;
1590 return *this;
1591 }
1592 BOOST_UBLAS_INLINE
1593 const_iterator1 &operator += (difference_type n) {
1594 it1_ += n;
1595 return *this;
1596 }
1597 BOOST_UBLAS_INLINE
1598 const_iterator1 &operator -= (difference_type n) {
1599 it1_ -= n;
1600 return *this;
1601 }
1602 BOOST_UBLAS_INLINE
1603 difference_type operator - (const const_iterator1 &it) const {
1604 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1605 return it1_ - it.it1_;
1606 }
1607
1608 // Dereference
1609 BOOST_UBLAS_INLINE
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 ())
1620 return *it1_;
1621#else
1622 size_type k = j;
1623 size_type l = (*this) ().upper () + i - j;
1624 if (k < (*this) ().size2 () &&
1625 l < (*this) ().lower () + 1 + (*this) ().upper ())
1626 return *it1_;
1627#endif
1628 return (*this) () (i, j);
1629 }
1630 BOOST_UBLAS_INLINE
1631 const_reference operator [] (difference_type n) const {
1632 return *(*this + n);
1633 }
1634
1635#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1636 BOOST_UBLAS_INLINE
1637#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1638 typename self_type::
1639#endif
1640 const_iterator2 begin () const {
1641 return (*this) ().find2 (1, index1 (), 0);
1642 }
1643 BOOST_UBLAS_INLINE
1644#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1645 typename self_type::
1646#endif
1647 const_iterator2 cbegin () const {
1648 return begin ();
1649 }
1650 BOOST_UBLAS_INLINE
1651#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1652 typename self_type::
1653#endif
1654 const_iterator2 end () const {
1655 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1656 }
1657 BOOST_UBLAS_INLINE
1658#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1659 typename self_type::
1660#endif
1661 const_iterator2 cend () const {
1662 return end ();
1663 }
1664 BOOST_UBLAS_INLINE
1665#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1666 typename self_type::
1667#endif
1668 const_reverse_iterator2 rbegin () const {
1669 return const_reverse_iterator2 (end ());
1670 }
1671 BOOST_UBLAS_INLINE
1672#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1673 typename self_type::
1674#endif
1675 const_reverse_iterator2 crbegin () const {
1676 return rbegin ();
1677 }
1678 BOOST_UBLAS_INLINE
1679#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1680 typename self_type::
1681#endif
1682 const_reverse_iterator2 rend () const {
1683 return const_reverse_iterator2 (begin ());
1684 }
1685 BOOST_UBLAS_INLINE
1686#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1687 typename self_type::
1688#endif
1689 const_reverse_iterator2 crend () const {
1690 return rend ();
1691 }
1692#endif
1693
1694 // Indices
1695 BOOST_UBLAS_INLINE
1696 size_type index1 () const {
1697 return it1_.index1 ();
1698 }
1699 BOOST_UBLAS_INLINE
1700 size_type index2 () const {
1701 return it1_.index2 ();
1702 }
1703
1704 // Assignment
1705 BOOST_UBLAS_INLINE
1706 const_iterator1 &operator = (const const_iterator1 &it) {
1707 container_const_reference<self_type>::assign (&it ());
1708 it1_ = it.it1_;
1709 return *this;
1710 }
1711
1712 // Comparison
1713 BOOST_UBLAS_INLINE
1714 bool operator == (const const_iterator1 &it) const {
1715 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1716 return it1_ == it.it1_;
1717 }
1718 BOOST_UBLAS_INLINE
1719 bool operator < (const const_iterator1 &it) const {
1720 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1721 return it1_ < it.it1_;
1722 }
1723
1724 private:
1725 const_subiterator1_type it1_;
1726 };
1727#endif
1728
1729 BOOST_UBLAS_INLINE
1730 const_iterator1 begin1 () const {
1731 return find1 (0, 0, 0);
1732 }
1733 BOOST_UBLAS_INLINE
1734 const_iterator1 cbegin1 () const {
1735 return begin1 ();
1736 }
1737 BOOST_UBLAS_INLINE
1738 const_iterator1 end1 () const {
1739 return find1 (0, size1 (), 0);
1740 }
1741 BOOST_UBLAS_INLINE
1742 const_iterator1 cend1 () const {
1743 return end1 ();
1744 }
1745
1746#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1747 class iterator1:
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> {
1752 public:
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;
1757
1758 typedef iterator2 dual_iterator_type;
1759 typedef reverse_iterator2 dual_reverse_iterator_type;
1760
1761 // Construction and destruction
1762 BOOST_UBLAS_INLINE
1763 iterator1 ():
1764 container_reference<self_type> (), it1_ () {}
1765 BOOST_UBLAS_INLINE
1766 iterator1 (self_type &m, const subiterator1_type &it1):
1767 container_reference<self_type> (m), it1_ (it1) {}
1768
1769 // Arithmetic
1770 BOOST_UBLAS_INLINE
1771 iterator1 &operator ++ () {
1772 ++ it1_;
1773 return *this;
1774 }
1775 BOOST_UBLAS_INLINE
1776 iterator1 &operator -- () {
1777 -- it1_;
1778 return *this;
1779 }
1780 BOOST_UBLAS_INLINE
1781 iterator1 &operator += (difference_type n) {
1782 it1_ += n;
1783 return *this;
1784 }
1785 BOOST_UBLAS_INLINE
1786 iterator1 &operator -= (difference_type n) {
1787 it1_ -= n;
1788 return *this;
1789 }
1790 BOOST_UBLAS_INLINE
1791 difference_type operator - (const iterator1 &it) const {
1792 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1793 return it1_ - it.it1_;
1794 }
1795
1796 // Dereference
1797 BOOST_UBLAS_INLINE
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 ())
1808 return *it1_;
1809#else
1810 size_type k = j;
1811 size_type l = (*this) ().upper () + i - j;
1812 if (k < (*this) ().size2 () &&
1813 l < (*this) ().lower () + 1 + (*this) ().upper ())
1814 return *it1_;
1815#endif
1816 return (*this) () (i, j);
1817 }
1818 BOOST_UBLAS_INLINE
1819 reference operator [] (difference_type n) const {
1820 return *(*this + n);
1821 }
1822
1823#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1824 BOOST_UBLAS_INLINE
1825#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1826 typename self_type::
1827#endif
1828 iterator2 begin () const {
1829 return (*this) ().find2 (1, index1 (), 0);
1830 }
1831 BOOST_UBLAS_INLINE
1832#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1833 typename self_type::
1834#endif
1835 iterator2 end () const {
1836 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1837 }
1838 BOOST_UBLAS_INLINE
1839#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1840 typename self_type::
1841#endif
1842 reverse_iterator2 rbegin () const {
1843 return reverse_iterator2 (end ());
1844 }
1845 BOOST_UBLAS_INLINE
1846#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1847 typename self_type::
1848#endif
1849 reverse_iterator2 rend () const {
1850 return reverse_iterator2 (begin ());
1851 }
1852#endif
1853
1854 // Indices
1855 BOOST_UBLAS_INLINE
1856 size_type index1 () const {
1857 return it1_.index1 ();
1858 }
1859 BOOST_UBLAS_INLINE
1860 size_type index2 () const {
1861 return it1_.index2 ();
1862 }
1863
1864 // Assignment
1865 BOOST_UBLAS_INLINE
1866 iterator1 &operator = (const iterator1 &it) {
1867 container_reference<self_type>::assign (&it ());
1868 it1_ = it.it1_;
1869 return *this;
1870 }
1871
1872 // Comparison
1873 BOOST_UBLAS_INLINE
1874 bool operator == (const iterator1 &it) const {
1875 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1876 return it1_ == it.it1_;
1877 }
1878 BOOST_UBLAS_INLINE
1879 bool operator < (const iterator1 &it) const {
1880 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1881 return it1_ < it.it1_;
1882 }
1883
1884 private:
1885 subiterator1_type it1_;
1886
1887 friend class const_iterator1;
1888 };
1889#endif
1890
1891 BOOST_UBLAS_INLINE
1892 iterator1 begin1 () {
1893 return find1 (0, 0, 0);
1894 }
1895 BOOST_UBLAS_INLINE
1896 iterator1 end1 () {
1897 return find1 (0, size1 (), 0);
1898 }
1899
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> {
1905 public:
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;
1912
1913 typedef const_iterator1 dual_iterator_type;
1914 typedef const_reverse_iterator1 dual_reverse_iterator_type;
1915
1916 // Construction and destruction
1917 BOOST_UBLAS_INLINE
1918 const_iterator2 ():
1919 container_const_reference<self_type> (), it2_ () {}
1920 BOOST_UBLAS_INLINE
1921 const_iterator2 (const self_type &m, const const_subiterator2_type &it2):
1922 container_const_reference<self_type> (m), it2_ (it2) {}
1923 BOOST_UBLAS_INLINE
1924 const_iterator2 (const iterator2 &it):
1925 container_const_reference<self_type> (it ()), it2_ (it.it2_) {}
1926
1927 // Arithmetic
1928 BOOST_UBLAS_INLINE
1929 const_iterator2 &operator ++ () {
1930 ++ it2_;
1931 return *this;
1932 }
1933 BOOST_UBLAS_INLINE
1934 const_iterator2 &operator -- () {
1935 -- it2_;
1936 return *this;
1937 }
1938 BOOST_UBLAS_INLINE
1939 const_iterator2 &operator += (difference_type n) {
1940 it2_ += n;
1941 return *this;
1942 }
1943 BOOST_UBLAS_INLINE
1944 const_iterator2 &operator -= (difference_type n) {
1945 it2_ -= n;
1946 return *this;
1947 }
1948 BOOST_UBLAS_INLINE
1949 difference_type operator - (const const_iterator2 &it) const {
1950 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1951 return it2_ - it.it2_;
1952 }
1953
1954 // Dereference
1955 BOOST_UBLAS_INLINE
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 ())
1966 return *it2_;
1967#else
1968 size_type k = j;
1969 size_type l = (*this) ().upper () + i - j;
1970 if (k < (*this) ().size2 () &&
1971 l < (*this) ().lower () + 1 + (*this) ().upper ())
1972 return *it2_;
1973#endif
1974 return (*this) () (i, j);
1975 }
1976 BOOST_UBLAS_INLINE
1977 const_reference operator [] (difference_type n) const {
1978 return *(*this + n);
1979 }
1980
1981#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1982 BOOST_UBLAS_INLINE
1983#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1984 typename self_type::
1985#endif
1986 const_iterator1 begin () const {
1987 return (*this) ().find1 (1, 0, index2 ());
1988 }
1989 BOOST_UBLAS_INLINE
1990#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1991 typename self_type::
1992#endif
1993 const_iterator1 cbegin () const {
1994 return begin ();
1995 }
1996 BOOST_UBLAS_INLINE
1997#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1998 typename self_type::
1999#endif
2000 const_iterator1 end () const {
2001 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2002 }
2003 BOOST_UBLAS_INLINE
2004#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2005 typename self_type::
2006#endif
2007 const_iterator1 cend () const {
2008 return end ();
2009 }
2010 BOOST_UBLAS_INLINE
2011#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2012 typename self_type::
2013#endif
2014 const_reverse_iterator1 rbegin () const {
2015 return const_reverse_iterator1 (end ());
2016 }
2017 BOOST_UBLAS_INLINE
2018#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2019 typename self_type::
2020#endif
2021 const_reverse_iterator1 crbegin () const {
2022 return rbegin ();
2023 }
2024 BOOST_UBLAS_INLINE
2025#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2026 typename self_type::
2027#endif
2028 const_reverse_iterator1 rend () const {
2029 return const_reverse_iterator1 (begin ());
2030 }
2031 BOOST_UBLAS_INLINE
2032#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2033 typename self_type::
2034#endif
2035 const_reverse_iterator1 crend () const {
2036 return rend ();
2037 }
2038#endif
2039
2040 // Indices
2041 BOOST_UBLAS_INLINE
2042 size_type index1 () const {
2043 return it2_.index1 ();
2044 }
2045 BOOST_UBLAS_INLINE
2046 size_type index2 () const {
2047 return it2_.index2 ();
2048 }
2049
2050 // Assignment
2051 BOOST_UBLAS_INLINE
2052 const_iterator2 &operator = (const const_iterator2 &it) {
2053 container_const_reference<self_type>::assign (&it ());
2054 it2_ = it.it2_;
2055 return *this;
2056 }
2057
2058 // Comparison
2059 BOOST_UBLAS_INLINE
2060 bool operator == (const const_iterator2 &it) const {
2061 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2062 return it2_ == it.it2_;
2063 }
2064 BOOST_UBLAS_INLINE
2065 bool operator < (const const_iterator2 &it) const {
2066 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2067 return it2_ < it.it2_;
2068 }
2069
2070 private:
2071 const_subiterator2_type it2_;
2072 };
2073#endif
2074
2075 BOOST_UBLAS_INLINE
2076 const_iterator2 begin2 () const {
2077 return find2 (0, 0, 0);
2078 }
2079 BOOST_UBLAS_INLINE
2080 const_iterator2 cbegin2 () const {
2081 return begin2 ();
2082 }
2083 BOOST_UBLAS_INLINE
2084 const_iterator2 end2 () const {
2085 return find2 (0, 0, size2 ());
2086 }
2087 BOOST_UBLAS_INLINE
2088 const_iterator2 cend2 () const {
2089 return end2 ();
2090 }
2091
2092#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2093 class iterator2:
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> {
2098 public:
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;
2103
2104 typedef iterator1 dual_iterator_type;
2105 typedef reverse_iterator1 dual_reverse_iterator_type;
2106
2107 // Construction and destruction
2108 BOOST_UBLAS_INLINE
2109 iterator2 ():
2110 container_reference<self_type> (), it2_ () {}
2111 BOOST_UBLAS_INLINE
2112 iterator2 (self_type &m, const subiterator2_type &it2):
2113 container_reference<self_type> (m), it2_ (it2) {}
2114
2115 // Arithmetic
2116 BOOST_UBLAS_INLINE
2117 iterator2 &operator ++ () {
2118 ++ it2_;
2119 return *this;
2120 }
2121 BOOST_UBLAS_INLINE
2122 iterator2 &operator -- () {
2123 -- it2_;
2124 return *this;
2125 }
2126 BOOST_UBLAS_INLINE
2127 iterator2 &operator += (difference_type n) {
2128 it2_ += n;
2129 return *this;
2130 }
2131 BOOST_UBLAS_INLINE
2132 iterator2 &operator -= (difference_type n) {
2133 it2_ -= n;
2134 return *this;
2135 }
2136 BOOST_UBLAS_INLINE
2137 difference_type operator - (const iterator2 &it) const {
2138 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2139 return it2_ - it.it2_;
2140 }
2141
2142 // Dereference
2143 BOOST_UBLAS_INLINE
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 ())
2154 return *it2_;
2155#else
2156 size_type k = j;
2157 size_type l = (*this) ().upper () + i - j;
2158 if (k < (*this) ().size2 () &&
2159 l < (*this) ().lower () + 1 + (*this) ().upper ())
2160 return *it2_;
2161#endif
2162 return (*this) () (i, j);
2163 }
2164 BOOST_UBLAS_INLINE
2165 reference operator [] (difference_type n) const {
2166 return *(*this + n);
2167 }
2168
2169#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2170 BOOST_UBLAS_INLINE
2171#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2172 typename self_type::
2173#endif
2174 iterator1 begin () const {
2175 return (*this) ().find1 (1, 0, index2 ());
2176 }
2177 BOOST_UBLAS_INLINE
2178#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2179 typename self_type::
2180#endif
2181 iterator1 end () const {
2182 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2183 }
2184 BOOST_UBLAS_INLINE
2185#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2186 typename self_type::
2187#endif
2188 reverse_iterator1 rbegin () const {
2189 return reverse_iterator1 (end ());
2190 }
2191 BOOST_UBLAS_INLINE
2192#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2193 typename self_type::
2194#endif
2195 reverse_iterator1 rend () const {
2196 return reverse_iterator1 (begin ());
2197 }
2198#endif
2199
2200 // Indices
2201 BOOST_UBLAS_INLINE
2202 size_type index1 () const {
2203 return it2_.index1 ();
2204 }
2205 BOOST_UBLAS_INLINE
2206 size_type index2 () const {
2207 return it2_.index2 ();
2208 }
2209
2210 // Assignment
2211 BOOST_UBLAS_INLINE
2212 iterator2 &operator = (const iterator2 &it) {
2213 container_reference<self_type>::assign (&it ());
2214 it2_ = it.it2_;
2215 return *this;
2216 }
2217
2218 // Comparison
2219 BOOST_UBLAS_INLINE
2220 bool operator == (const iterator2 &it) const {
2221 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2222 return it2_ == it.it2_;
2223 }
2224 BOOST_UBLAS_INLINE
2225 bool operator < (const iterator2 &it) const {
2226 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2227 return it2_ < it.it2_;
2228 }
2229
2230 private:
2231 subiterator2_type it2_;
2232
2233 friend class const_iterator2;
2234 };
2235#endif
2236
2237 BOOST_UBLAS_INLINE
2238 iterator2 begin2 () {
2239 return find2 (0, 0, 0);
2240 }
2241 BOOST_UBLAS_INLINE
2242 iterator2 end2 () {
2243 return find2 (0, 0, size2 ());
2244 }
2245
2246 // Reverse iterators
2247
2248 BOOST_UBLAS_INLINE
2249 const_reverse_iterator1 rbegin1 () const {
2250 return const_reverse_iterator1 (end1 ());
2251 }
2252 BOOST_UBLAS_INLINE
2253 const_reverse_iterator1 crbegin1 () const {
2254 return rbegin1 ();
2255 }
2256 BOOST_UBLAS_INLINE
2257 const_reverse_iterator1 rend1 () const {
2258 return const_reverse_iterator1 (begin1 ());
2259 }
2260 BOOST_UBLAS_INLINE
2261 const_reverse_iterator1 crend1 () const {
2262 return rend1 ();
2263 }
2264
2265 BOOST_UBLAS_INLINE
2266 reverse_iterator1 rbegin1 () {
2267 return reverse_iterator1 (end1 ());
2268 }
2269 BOOST_UBLAS_INLINE
2270 reverse_iterator1 rend1 () {
2271 return reverse_iterator1 (begin1 ());
2272 }
2273
2274 BOOST_UBLAS_INLINE
2275 const_reverse_iterator2 rbegin2 () const {
2276 return const_reverse_iterator2 (end2 ());
2277 }
2278 BOOST_UBLAS_INLINE
2279 const_reverse_iterator2 crbegin2 () const {
2280 return rbegin2 ();
2281 }
2282 BOOST_UBLAS_INLINE
2283 const_reverse_iterator2 rend2 () const {
2284 return const_reverse_iterator2 (begin2 ());
2285 }
2286 BOOST_UBLAS_INLINE
2287 const_reverse_iterator2 crend2 () const {
2288 return rend2 ();
2289 }
2290
2291 BOOST_UBLAS_INLINE
2292 reverse_iterator2 rbegin2 () {
2293 return reverse_iterator2 (end2 ());
2294 }
2295 BOOST_UBLAS_INLINE
2296 reverse_iterator2 rend2 () {
2297 return reverse_iterator2 (begin2 ());
2298 }
2299
2300 private:
2301 matrix_closure_type data_;
2302 size_type lower_;
2303 size_type upper_;
2304 typedef const value_type const_value_type;
2305 static const_value_type zero_;
2306 };
2307
2308 // Specialization for temporary_traits
2309 template <class M>
2310 struct vector_temporary_traits< banded_adaptor<M> >
2311 : vector_temporary_traits< M > {} ;
2312 template <class M>
2313 struct vector_temporary_traits< const banded_adaptor<M> >
2314 : vector_temporary_traits< M > {} ;
2315
2316 template <class M>
2317 struct matrix_temporary_traits< banded_adaptor<M> >
2318 : matrix_temporary_traits< M > {} ;
2319 template <class M>
2320 struct matrix_temporary_traits< const banded_adaptor<M> >
2321 : matrix_temporary_traits< M > {} ;
2322
2323
2324 template<class M>
2325 typename banded_adaptor<M>::const_value_type banded_adaptor<M>::zero_ = value_type/*zero*/();
2326
2327 /** \brief A diagonal matrix adaptator: convert a any matrix into a diagonal matrix expression
2328 *
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$.
2331 *
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.
2335 *
2336 * \tparam M the type of matrix used to generate the diagonal matrix
2337 */
2338
2339 template<class M>
2340 class diagonal_adaptor:
2341 public banded_adaptor<M> {
2342 public:
2343 typedef M matrix_type;
2344 typedef banded_adaptor<M> adaptor_type;
2345
2346 // Construction and destruction
2347 BOOST_UBLAS_INLINE
2348 diagonal_adaptor ():
2349 adaptor_type () {}
2350 BOOST_UBLAS_INLINE
2351 diagonal_adaptor (matrix_type &data):
2352 adaptor_type (data) {}
2353 BOOST_UBLAS_INLINE
2354 ~diagonal_adaptor () {}
2355
2356 // Assignment
2357 BOOST_UBLAS_INLINE
2358 diagonal_adaptor &operator = (const diagonal_adaptor &m) {
2359 adaptor_type::operator = (m);
2360 return *this;
2361 }
2362 template<class AE>
2363 BOOST_UBLAS_INLINE
2364 diagonal_adaptor &operator = (const matrix_expression<AE> &ae) {
2365 adaptor_type::operator = (ae);
2366 return *this;
2367 }
2368 };
2369
2370}}}
2371
2372#endif