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