00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #ifndef SCPPNT_REGION2D_H
00044 #define SCPPNT_REGION2D_H
00045
00046 #ifdef SCPPNT_NO_DIR_PREFIX
00047 #include "scppnt.h"
00048 #include "slice_pointer.h"
00049 #include "scppnt_error.h"
00050 #include "matop.h"
00051 #include "index.h"
00052 #else
00053 #include "scppnt/scppnt.h"
00054 #include "scppnt/slice_pointer.h"
00055 #include "scppnt/scppnt_error.h"
00056 #include "scppnt/matop.h"
00057 #include "scppnt/index.h"
00058 #endif
00059
00060 #ifndef SCPPNT_NO_IO
00061 #include <iostream>
00062 #endif
00063
00064 namespace SCPPNT
00065 {
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095 template<class T, class IT2D> class Region2D_iterator
00096 {
00097 public:
00098
00099
00100 typedef T value_type;
00101 typedef std::forward_iterator_tag iterator_category;
00102 typedef Subscript difference_type;
00103
00104 typedef IT2D iterator_2D;
00105
00106
00107 typedef typename std::iterator_traits<iterator_2D>::value_type iterator_1D;
00108
00109 typedef typename std::iterator_traits<iterator_1D>::pointer pointer;
00110
00111 typedef typename std::iterator_traits<iterator_1D>::reference reference;
00112
00113
00114 Region2D_iterator(IT2D it, Subscript n2d, Subscript n1d, Subscript initial2 = 1,
00115 Subscript initial1 = 1);
00116
00117
00118 reference operator*() const
00119 {
00120 return *current_;
00121 }
00122
00123
00124 Region2D_iterator& operator++()
00125 {
00126 ++current_;
00127 if (current_ == end_)
00128 {
00129 ++i2D_;
00130 if (i2D_ != end2D_)
00131 {
00132 current_ = *i2D_;
00133 end_ = current_ + n1D_;
00134 }
00135 }
00136 return *this;
00137 }
00138
00139
00140 Region2D_iterator operator++(int)
00141 {
00142 Region2D_iterator t(*this);
00143 ++current_;
00144 if (current_ == end_)
00145 {
00146 ++i2D_;
00147 if (i2D_ != end2D_)
00148 {
00149 current_ = *i2D_;
00150 end_ = current_ + n1D_;
00151 }
00152 }
00153 return t;
00154 }
00155
00156 #ifdef SCPPNT_MEMBER_COMPARISONS
00157
00158
00159 bool operator==(const Region2D_iterator<T, IT2D> &rhs)
00160 {
00161 return current_ == rhs.current_;
00162 }
00163
00164
00165 bool operator!=(const Region2D_iterator<T, IT2D> &rhs)
00166 {
00167 return current_ != rhs.current_;
00168 }
00169
00170 #else
00171
00172
00173 friend bool operator==(const Region2D_iterator<T, IT2D> &lhs,
00174 const Region2D_iterator<T, IT2D> &rhs)
00175 {
00176 return lhs.current_ == rhs.current_;
00177 }
00178
00179
00180 friend bool operator!=(const Region2D_iterator<T, IT2D> &lhs,
00181 const Region2D_iterator<T, IT2D> &rhs)
00182 {
00183 return lhs.current_ != rhs.current_;
00184 }
00185
00186 #endif
00187
00188 private:
00189
00190
00191 iterator_2D i2D_;
00192
00193
00194 iterator_1D current_;
00195
00196
00197 iterator_1D end_;
00198
00199
00200 iterator_2D end2D_;
00201
00202
00203 Subscript n1D_;
00204
00205 };
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 template<class T, class IT2D> Region2D_iterator<T, IT2D>::Region2D_iterator(IT2D it,
00217 Subscript n2d, Subscript n1d, Subscript initial2, Subscript initial1) :
00218 i2D_(it), n1D_(n1d)
00219 {
00220 if (initial2 >= n2d || initial1 >= n1d)
00221 throw BadDimension("SCPPNT::Region2D_iterator::region2D_iterator");
00222 end2D_ = it + n2d;
00223 i2D_ += initial2 - 1;
00224 end_ = *i2D_ + n1D_;
00225 current_ = *i2D_ + initial1 - 1;
00226 }
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241 template<class Array2D> class Region2D
00242 {
00243 public:
00244
00245 typedef Array2D array_type;
00246 typedef typename Array2D::size_type size_type;
00247 typedef typename Array2D::value_type value_type;
00248 typedef typename Array2D::element_type element_type;
00249 typedef typename Array2D::pointer pointer;
00250 typedef typename Array2D::reference reference;
00251 typedef typename Array2D::const_reference const_reference;
00252
00253
00254
00255
00256 typedef typename Array2D::row_iterator row_iterator;
00257
00258
00259 typedef typename Array2D::const_row_iterator const_row_iterator;
00260
00261
00262 typedef typename Array2D::column_iterator column_iterator;
00263
00264
00265 typedef typename Array2D::const_column_iterator const_column_iterator;
00266
00267
00268
00269
00270
00271 #ifdef SCPPNT_BOUNDS_CHECK
00272 typedef slice_pointer_base<row_iterator, row_iterator*, row_iterator&> rows_iterator;
00273 typedef slice_pointer_base<column_iterator, column_iterator*, column_iterator&>
00274 columns_iterator;
00275 typedef
00276 slice_pointer_base<const_row_iterator, const_row_iterator*, const_row_iterator&>
00277 const_rows_iterator;
00278 typedef
00279 slice_pointer_base<const_column_iterator, const_column_iterator*, const_column_iterator&>
00280 const_columns_iterator;
00281 #else
00282
00283 typedef row_iterator* rows_iterator;
00284
00285
00286 typedef const_row_iterator* const_rows_iterator;
00287
00288
00289 typedef column_iterator* columns_iterator;
00290
00291
00292 typedef const_column_iterator* const_columns_iterator;
00293 #endif
00294
00295
00296 typedef Region2D_iterator<value_type, typename Region2D<Array2D>::rows_iterator> iterator;
00297
00298 typedef Region2D_iterator<value_type, typename Region2D<Array2D>::const_rows_iterator>
00299 const_iterator;
00300
00301
00302 typedef typename Array2D::diag_iterator diag_iterator;
00303
00304 typedef typename Array2D::const_diag_iterator const_diag_iterator;
00305
00306
00307 const array_type & array() const
00308 {
00309 return A_;
00310 }
00311
00312
00313 Subscript num_rows() const
00314 {
00315 return dim_[0];
00316 }
00317
00318
00319 Subscript num_columns() const
00320 {
00321 return dim_[1];
00322 }
00323
00324
00325 Subscript lbound() const
00326 {
00327 return A_.lbound();
00328 }
00329
00330
00331 Subscript dim(Subscript i) const
00332 {
00333 if (i< 1 || i> 2)throw InvalidArgument("argument must be 1 or 2", "SCPPNT::Region2D::dim");
00334 return (i==1) ? dim_[0] : dim_[1];
00335 }
00336
00337
00338 Subscript size() const
00339 {
00340 return dim_[0]*dim_[1];
00341 }
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 Region2D(Array2D &A, Subscript i1, Subscript i2, Subscript j1, Subscript j2) : A_(A)
00355 {
00356 initialize(i1, i2, j1, j2);
00357 }
00358
00359
00360
00361
00362
00363
00364
00365 Region2D(Array2D &A, const Index1D &I, const Index1D &J) : A_(A)
00366 {
00367 initialize(I.lbound(), I.ubound(), J.lbound(), J.ubound());
00368 }
00369
00370
00371 Region2D(Array2D &A) : A_(A)
00372 {
00373 initialize(1, A.num_rows(), 1, A.num_columns());
00374 }
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384 Region2D(Region2D<Array2D> &A, Subscript i1, Subscript i2, Subscript j1, Subscript j2)
00385 : A_(A.A_)
00386 {
00387 initialize(i1 + A.offset_[0], i2 + A.offset_[0], j1 + A.offset_[1], j2 + A.offset_[1]);
00388 }
00389
00390
00391 Region2D(const Region2D<Array2D> &A) : A_(A.A_)
00392 {
00393 initialize(1 + A.offset_[0], A.dim_[0] + A.offset_[0],
00394 1 + A.offset_[1], A.dim_[1] + A.offset_[1]);
00395 }
00396
00397
00398 ~Region2D()
00399 {
00400 if (row_ != 0) delete [] (row_);
00401 if (column_ != 0) delete [] (column_);
00402 if (const_row_ != 0) delete [] (const_row_);
00403 if (const_column_ != 0) delete [] (const_column_);}
00404
00405
00406
00407
00408 row_iterator operator[](Subscript i)
00409 {
00410 return row_[i];
00411 }
00412
00413
00414 const_row_iterator operator[](Subscript i) const
00415 {
00416 return row_[i];
00417 }
00418
00419
00420 reference operator()(Subscript i, Subscript j)
00421 {
00422 return row_[i-1][j-1];
00423 }
00424
00425
00426 const_reference operator()(Subscript i, Subscript j) const
00427 {
00428 return row_[i-1][j-1];
00429 }
00430
00431
00432
00433
00434
00435
00436
00437
00438 Region2D<Array2D> operator()(Subscript i1, Subscript i2,
00439 Subscript j1, Subscript j2)
00440 {
00441 return Region2D<Array2D>(A_,
00442 i1+offset_[0], offset_[0] + i2,
00443 j1+offset_[1], offset_[1] + j2);
00444 }
00445
00446
00447
00448
00449
00450
00451 Region2D<Array2D> operator()(const Index1D &I,
00452 const Index1D &J)
00453 {
00454 return Region2D<Array2D>(A_, I.lbound()+offset_[0],
00455 offset_[0] + I.ubound(), offset_[1]+J.lbound(),
00456 offset_[1] + J.ubound());
00457 }
00458
00459
00460
00461
00462 rows_iterator begin_rows();
00463
00464
00465 const_rows_iterator begin_rows() const;
00466
00467
00468 columns_iterator begin_columns();
00469
00470
00471 const_columns_iterator begin_columns() const;
00472
00473
00474 rows_iterator end_rows();
00475
00476
00477 const_rows_iterator end_rows() const;
00478
00479
00480 columns_iterator end_columns();
00481
00482
00483 const_columns_iterator end_columns() const;
00484
00485
00486
00487
00488
00489 row_iterator begin_row(Subscript index)
00490 {
00491 return row_[index-1];
00492 }
00493
00494
00495 const_row_iterator begin_row(Subscript index) const
00496 {
00497 return const_row_[index-1];
00498 }
00499
00500
00501 column_iterator begin_column(Subscript index)
00502 {
00503 return column_[index-1];
00504 }
00505
00506
00507 const_column_iterator begin_column(Subscript index) const
00508 {
00509 return const_column_[index-1];
00510 }
00511
00512
00513 row_iterator end_row(Subscript index)
00514 {
00515 return row_[index];
00516 }
00517
00518
00519 const_row_iterator end_row(Subscript index) const
00520 {
00521 return const_row_[index];
00522 }
00523
00524
00525 column_iterator end_column(Subscript index)
00526 {
00527 return column_[index];
00528 }
00529
00530
00531 const_column_iterator end_column(Subscript index) const
00532 {
00533 return const_column_[index];
00534 }
00535
00536
00537
00538
00539
00540
00541 iterator begin()
00542 {
00543 return iterator(begin_rows(), num_rows(), num_columns());
00544 }
00545
00546
00547 const_iterator begin() const
00548 {
00549 return const_iterator(begin_rows(), num_rows(), num_columns());
00550 }
00551
00552
00553 iterator end()
00554 {
00555 return iterator(begin_rows(), num_rows(), num_columns(), dim_[0]+1, 1);
00556 }
00557
00558
00559 const_iterator end() const
00560 {
00561 return const_iterator(begin_rows(), num_rows(), num_columns(), dim_[0]+1, 1);
00562 }
00563
00564
00565
00566
00567
00568 diag_iterator begin_diagonal(Subscript row, Subscript column);
00569
00570
00571
00572 const_diag_iterator begin_diagonal(Subscript row, Subscript column) const;
00573
00574
00575
00576 diag_iterator end_diagonal(Subscript row, Subscript column);
00577
00578
00579
00580 const_diag_iterator end_diagonal(Subscript row, Subscript column) const;
00581
00582
00583
00584
00585
00586 Region2D<Array2D>& operator=(const value_type& scalar)
00587 {
00588 iterator it = begin();
00589 for (Subscript j = size(); j--; ++it)
00590 *it = scalar;
00591 }
00592
00593
00594 Region2D<Array2D> & operator=(const Region2D<Array2D> &R);
00595
00596
00597
00598
00599
00600
00601 template<class MAT>
00602 Region2D<Array2D>& operator+=(const MAT &rhs)
00603 {
00604 matadd(*this, rhs);
00605 return *this;
00606 }
00607
00608
00609 template<class MAT>
00610 Region2D<Array2D>& operator-=(const MAT &rhs)
00611 {
00612 matsub(*this, rhs);
00613 return *this;
00614 }
00615
00616
00617 template<class MAT>
00618 Region2D<Array2D>& operator*=(const MAT &rhs)
00619 {
00620 matmult_assign(*this, rhs);
00621 return *this;
00622 }
00623
00624
00625 Region2D<Array2D>& operator+=(const value_type& value)
00626 { iterator it = begin();
00627 for (Subscript j = size(); j--; ++it)
00628 *it += value;
00629 }
00630
00631
00632 Region2D<Array2D>& operator-=(const value_type& value)
00633 {
00634 iterator it = begin();
00635 for (Subscript j = size(); j--; ++it)
00636 *it -= value;
00637 }
00638
00639
00640 Region2D<Array2D>& operator*=(const value_type& value)
00641 { iterator it = begin(); for (Subscript j = size(); j--; ++it) *it *= value;}
00642
00643
00644 Region2D<Array2D>& operator/=(const value_type& value)
00645 {
00646 iterator it = begin();
00647 for (Subscript j = size(); j--; ++it)
00648 *it /= value;
00649 }
00650
00651 private:
00652
00653
00654 void initialize(Subscript rowLow, Subscript rowHi, Subscript colLow, Subscript colHi);
00655
00656 Array2D& A_;
00657 Subscript offset_[2];
00658 Subscript dim_[2];
00659
00660 row_iterator *row_;
00661 const_row_iterator *const_row_;
00662
00663
00664 column_iterator *column_;
00665 const_column_iterator *const_column_;
00666
00667
00668
00669 Subscript diagonal_size(Subscript row, Subscript column) const;
00670
00671 };
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681 template <class Array2D>
00682 void Region2D<Array2D>::initialize(Subscript rowLow, Subscript rowHi,
00683 Subscript colLow, Subscript colHi)
00684 {
00685 Subscript lowerBound = A_.lbound();
00686
00687 if ( rowHi < rowLow || rowLow < lowerBound || rowHi > A_.num_rows() )
00688 throw RuntimeError("Bad row index range", "SCPPNT::Region2D<Array2D>::initialize()");
00689 if ( colHi < colLow || colLow < lowerBound || colHi > A_.num_columns() )
00690 throw RuntimeError("Bad column index range", "SCPPNT::Region2D<Array2D>::initialize()");
00691
00692 offset_[0] = rowLow-lowerBound;
00693 offset_[1] = colLow-lowerBound;
00694 dim_[0] = rowHi-rowLow+1;
00695 dim_[1] = colHi-colLow+1;
00696
00697
00698 row_ = new row_iterator[dim_[0]+1];
00699 const_row_ = new const_row_iterator[dim_[0]+1];
00700
00701
00702 Subscript i;
00703 for (i=0; i<dim_[0]; i++)
00704 {
00705 row_[i] = A_.begin_row(i+offset_[0]+1) + offset_[1];
00706 #ifdef SCPPNT_BOUNDS_CHECK
00707 row_[i].set_size(dim_[1]);
00708 #endif
00709 const_row_[i] = row_[i];
00710 }
00711
00712
00713 row_[i] = row_[i-1] + offset_[1];
00714 const_row_[i] = row_[i];
00715 #ifdef SCPPNT_BOUNDS_CHECK
00716 row_[i].set_size(1);
00717 #endif
00718
00719
00720 column_ = new column_iterator[dim_[1]+1];
00721 const_column_ = new const_column_iterator[dim_[1]+1];
00722
00723
00724 for (i=0; i<dim_[1]; ++i)
00725 {
00726 column_[i] = A_.begin_column(i+offset_[1]+1) + offset_[0];
00727 #ifdef SCPPNT_BOUNDS_CHECK
00728 column_[i].set_size(dim_[0]);
00729 #endif
00730 const_column_[i] = column_[i];
00731 }
00732
00733
00734 column_[i] = column_[i-1] + offset_[0];
00735 const_column_[i] = column_[i];
00736 #ifdef SCPPNT_BOUNDS_CHECK
00737 column_[i].set_size(1);
00738 #endif
00739
00740 }
00741
00742
00743 template <class Array2D>
00744 inline typename Region2D<Array2D>::rows_iterator Region2D<Array2D>::begin_rows()
00745 {
00746 #ifdef SCPPNT_BOUNDS_CHECK
00747 return rows_iterator(row_, 1, dim_[0]);
00748 #else
00749 return row_;
00750 #endif
00751 }
00752
00753
00754 template <class Array2D>
00755 inline typename Region2D<Array2D>::const_rows_iterator Region2D<Array2D>::begin_rows()
00756 const
00757 {
00758 #ifdef SCPPNT_BOUNDS_CHECK
00759 return const_rows_iterator(const_row_, 1, dim_[0]);
00760 #else
00761 return const_row_;
00762 #endif
00763 }
00764
00765
00766 template <class Array2D>
00767 inline typename Region2D<Array2D>::columns_iterator Region2D<Array2D>::begin_columns()
00768 {
00769 #ifdef SCPPNT_BOUNDS_CHECK
00770 return columns_iterator(column_, 1, dim_[1]);
00771 #else
00772 return column_;
00773 #endif
00774 }
00775
00776
00777 template <class Array2D>
00778 inline typename Region2D<Array2D>::const_columns_iterator Region2D<Array2D>::begin_columns() const
00779 {
00780 #ifdef SCPPNT_BOUNDS_CHECK
00781 return const_columns_iterator(const_column_, 1, dim_[1]);
00782 #else
00783 return const_column_;
00784 #endif
00785 }
00786
00787
00788 template <class Array2D>
00789 inline typename Region2D<Array2D>::rows_iterator Region2D<Array2D>::end_rows()
00790 {
00791 #ifdef SCPPNT_BOUNDS_CHECK
00792 return rows_iterator(row_ + dim_[0], 1, 1);
00793 #else
00794 return row_ + dim_[0];
00795 #endif
00796 }
00797
00798
00799 template <class Array2D>
00800 inline typename Region2D<Array2D>::const_rows_iterator Region2D<Array2D>::end_rows() const
00801 {
00802 #ifdef SCPPNT_BOUNDS_CHECK
00803 return const_rows_iterator(const_row_ + dim_[0], 1, 1);
00804 #else
00805 return (const_rows_iterator) row_ + dim_[0];
00806 #endif
00807 }
00808
00809
00810 template <class Array2D>
00811 inline typename Region2D<Array2D>::columns_iterator Region2D<Array2D>::end_columns()
00812 {
00813 #ifdef SCPPNT_BOUNDS_CHECK
00814 return columns_iterator(column_ + dim_[1], 1, 1);
00815 #else
00816 return column_ + dim_[1];
00817 #endif
00818 }
00819
00820
00821 template <class Array2D>
00822 inline typename Region2D<Array2D>::const_columns_iterator Region2D<Array2D>::end_columns() const
00823 {
00824 #ifdef SCPPNT_BOUNDS_CHECK
00825 return const_columns_iterator(const_column_ + dim_[1], 1, 1);
00826 #else
00827 return const_column_ + dim_[1];
00828 #endif
00829 }
00830
00831
00832 template <class Array2D>
00833 inline Subscript Region2D<Array2D>::diagonal_size(Subscript row, Subscript column) const
00834 {
00835 Subscript ncolumn = dim_[0] - column + 1;
00836 Subscript nrow = dim_[1] - row + 1;
00837 return (ncolumn > nrow) ? nrow : ncolumn;
00838
00839 }
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852 template <class Array2D>
00853 inline typename Region2D<Array2D>::diag_iterator Region2D<Array2D>::begin_diagonal(Subscript row, Subscript column)
00854 {
00855
00856 diag_iterator it = A_.begin_diagonal(row + offset_[0], column + offset_[1]);
00857
00858 #ifdef SCPPNT_BOUNDS_CHECK
00859 it.set_size(diagonal_size(row, column));
00860 #endif
00861
00862 return it;
00863 }
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875 template <class Array2D>
00876 inline typename Region2D<Array2D>::const_diag_iterator Region2D<Array2D>::begin_diagonal(Subscript row, Subscript column) const
00877 {
00878 const_diag_iterator it = A_.begin_diagonal(row + offset_[0], column + offset_[1]);
00879
00880 #ifdef SCPPNT_BOUNDS_CHECK
00881 it.set_size(diagonal_size(row, column));
00882 #endif
00883
00884 return it;
00885 }
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898 template <class Array2D>
00899 inline typename Region2D<Array2D>::diag_iterator Region2D<Array2D>::end_diagonal(Subscript row, Subscript column)
00900 {
00901 diag_iterator it = begin_diagonal(row, column);
00902 return it + diagonal_size(row, column);
00903 }
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916 template <class Array2D>
00917 inline typename Region2D<Array2D>::const_diag_iterator Region2D<Array2D>::end_diagonal(Subscript row, Subscript column) const
00918 {
00919 const_diag_iterator it = begin_diagonal(row, column);
00920 return it + diagonal_size(row, column);
00921 }
00922
00923
00924
00925
00926
00927
00928
00929 template <class Array2D>
00930 Region2D<Array2D> & Region2D<Array2D>::operator=(const Region2D<Array2D> &R)
00931 {
00932 Subscript M = num_rows();
00933 Subscript N = num_columns();
00934
00935 if (M != R.num_rows() || N != R.num_columns())
00936 {
00937 throw BadDimension("Region2D<Array2D>::operator=(const Region2D)");
00938 }
00939
00940 const_iterator ir = R.begin();
00941 iterator ithis = begin();
00942 for (int i = M*N; i--; ++ir, ++ithis)
00943 {
00944 *ithis = *ir;
00945 }
00946 return *this;
00947 }
00948
00949 #ifndef SCPPNT_NO_IO
00950
00951
00952
00953
00954
00955
00956
00957 template <class Array2D>
00958 std::ostream& operator<<(std::ostream &s, const Region2D<Array2D> &A)
00959 {
00960 Subscript M=A.num_rows();
00961 Subscript N=A.num_columns();
00962
00963 Subscript start = A.lbound();
00964 Subscript Mend = M + A.lbound() - 1;
00965 Subscript Nend = N + A.lbound() - 1;
00966
00967 s << M << " " << N << std::endl;
00968 for (Subscript i=start; i<=Mend; i++)
00969 {
00970 for (Subscript j=start; j<=Nend; j++)
00971 {
00972 s << A(i,j) << " ";
00973 }
00974 s << std::endl;
00975 }
00976 return s;
00977 }
00978 #endif
00979
00980
00981
00982
00983 template <class Array2D, class M>
00984 inline M operator+(const Region2D<Array2D> &A,
00985 const M &B)
00986 {
00987 M sum(A.num_rows(), A.num_columns());
00988 matadd(A, B, sum);
00989 return sum;
00990 }
00991
00992
00993 template <class Array2D, class M>
00994 inline M operator-(const Region2D<Array2D> &A,
00995 const M &B)
00996 {
00997 M diff(A.num_rows(), A.num_columns());
00998 matsub(A, B, diff);
00999 return diff;
01000 }
01001
01002
01003 template <class Array2D, class M>
01004 inline M operator*(const Region2D<Array2D> &A,
01005 const M &B)
01006 {
01007 M c(A.num_rows(), B.num_columns());
01008 matmult(c, A, B);
01009 return c;
01010 }
01011
01012
01013 template <class T, class Array2D>
01014 inline Vector<T> operator*(const Region2D<Array2D> &A, const Vector<T> &x)
01015 {
01016 return matrix_times_vector< Region2D<Array2D>, Vector<T> >(A, x);
01017 }
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032 template <class Array2D>
01033 class const_Region2D
01034 {
01035 public:
01036
01037 typedef const Array2D array_type;
01038 typedef typename Array2D::size_type size_type;
01039 typedef typename Array2D::value_type value_type;
01040 typedef typename Array2D::element_type element_type;
01041 typedef typename Array2D::pointer pointer;
01042 typedef typename Array2D::reference reference;
01043 typedef typename Array2D::const_reference const_reference;
01044
01045
01046
01047
01048 typedef typename Array2D::const_row_iterator const_row_iterator;
01049
01050 typedef typename Array2D::const_column_iterator const_column_iterator;
01051
01052
01053
01054
01055
01056
01057 #ifdef SCPPNT_BOUNDS_CHECK
01058 typedef slice_pointer_base<const_row_iterator, const_row_iterator*, const_row_iterator&>
01059 const_rows_iterator;
01060 typedef slice_pointer_base<const_column_iterator, const_column_iterator*, const_column_iterator&>
01061 const_columns_iterator;
01062 #else
01063
01064 typedef const_row_iterator* const_rows_iterator;
01065
01066
01067 typedef const_column_iterator* const_columns_iterator;
01068 #endif
01069
01070
01071 typedef Region2D_iterator<value_type, typename const_Region2D<Array2D>::const_rows_iterator>
01072 const_iterator;
01073
01074
01075 typedef typename Array2D::const_diag_iterator const_diag_iterator;
01076
01077
01078 const array_type & array() const
01079 {
01080 return A_;
01081 }
01082
01083
01084 Subscript num_rows() const
01085 {
01086 return dim_[0];
01087 }
01088
01089
01090 Subscript num_columns() const
01091 {
01092 return dim_[1];
01093 }
01094
01095
01096 Subscript lbound() const
01097 {
01098 return A_.lbound();
01099 }
01100
01101
01102 Subscript dim(Subscript i) const
01103 {
01104 if (i < 1 || i > 2) throw InvalidArgument("argument must be 1 or 2", "SCPPNT::const_Region2D::dim");
01105 return (i==1) ? dim_[0] : dim_[1];
01106 }
01107
01108
01109 Subscript size() const
01110 {
01111 return dim_[0]*dim_[1];
01112 }
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125 const_Region2D(const Array2D &A, Subscript i1, Subscript i2, Subscript j1,
01126 Subscript j2) : A_(A)
01127 {
01128 initialize(i1, i2, j1, j2);
01129 }
01130
01131
01132
01133
01134
01135
01136
01137 const_Region2D(const Array2D &A, const Index1D &I, const Index1D &J) : A_(A)
01138 {
01139 initialize(I.lbound(), I.ubound(), J.lbound(), J.ubound());
01140 }
01141
01142
01143 const_Region2D(const Array2D &A) : A_(A)
01144 {
01145 initialize(1, A.num_rows(), 1, A.num_columns());
01146 }
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156 const_Region2D(const const_Region2D<Array2D> &A, Subscript i1, Subscript i2, Subscript j1, Subscript j2)
01157 : A_(A.A_)
01158 {
01159 initialize(i1 + A.offset_[0], i2 + A.offset_[0], j1 + A.offset_[1], j2 + A.offset_[1]);
01160 }
01161
01162
01163 const_Region2D(const const_Region2D<Array2D> &A) : A_(A.A_)
01164 {
01165 initialize(1 + A.offset_[0], A.dim_[0] + A.offset_[0], 1 + A.offset_[1], A.dim_[1] + A.offset_[1]);
01166 }
01167
01168
01169 ~const_Region2D()
01170 {
01171 if (const_row_ != 0) delete [] (const_row_);
01172 if (const_column_ != 0) delete [] (const_column_);
01173 }
01174
01175
01176
01177
01178 const_row_iterator operator[](Subscript i) const
01179 {
01180 return const_row_[i];
01181 }
01182
01183
01184 const_reference operator()(Subscript i, Subscript j) const
01185 {
01186 return const_row_[i-1][j-1];
01187 }
01188
01189
01190
01191
01192
01193
01194
01195
01196 const_Region2D<Array2D> operator()(Subscript i1, Subscript i2,
01197 Subscript j1, Subscript j2)
01198 {
01199 return const_Region2D<Array2D>(A_,
01200 i1+offset_[0], offset_[0] + i2,
01201 j1+offset_[1], offset_[1] + j2);
01202 }
01203
01204
01205
01206
01207
01208
01209 const_Region2D<Array2D> operator()(const Index1D &I,
01210 const Index1D &J)
01211 {
01212 return const_Region2D<Array2D>(A_, I.lbound()+offset_[0],
01213 offset_[0] + I.ubound(), offset_[1]+J.lbound(),
01214 offset_[1] + J.ubound());
01215 }
01216
01217
01218
01219
01220 const_rows_iterator begin_rows() const;
01221
01222
01223 const_columns_iterator begin_columns() const;
01224
01225
01226 const_rows_iterator end_rows() const;
01227
01228
01229 const_columns_iterator end_columns() const;
01230
01231
01232
01233
01234
01235
01236
01237
01238 const_row_iterator begin_row(Subscript index) const
01239 {
01240 return const_row_[index-1];
01241 }
01242
01243
01244 const_column_iterator begin_column(Subscript index) const
01245 {
01246 return const_column_[index-1];
01247 }
01248
01249
01250 const_row_iterator end_row(Subscript index) const
01251 {
01252 return const_row_[index];
01253 }
01254
01255
01256 const_column_iterator end_column(Subscript index) const
01257 {
01258 return const_column_[index];
01259 }
01260
01261
01262
01263
01264
01265
01266 const_iterator begin() const
01267 {
01268 return const_iterator(begin_rows(), num_rows(), num_columns());
01269 }
01270
01271
01272 const_iterator end() const
01273 {
01274 return const_iterator(begin_rows(), num_rows(), num_columns(), dim_[0]+1, 1);
01275 }
01276
01277
01278
01279
01280 const_diag_iterator begin_diagonal(Subscript row, Subscript column) const;
01281
01282
01283 const_diag_iterator end_diagonal(Subscript row, Subscript column) const;
01284
01285 private:
01286
01287
01288 void initialize(Subscript rowLow, Subscript rowHi, Subscript colLow, Subscript colHi);
01289
01290 const Array2D& A_;
01291 Subscript offset_[2];
01292 Subscript dim_[2];
01293
01294 const_row_iterator *const_row_;
01295
01296 const_column_iterator *const_column_;
01297
01298
01299 Subscript diagonal_size(Subscript row, Subscript column) const;
01300 };
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311 template <class Array2D>
01312 void const_Region2D<Array2D>::initialize(Subscript rowLow, Subscript rowHi,
01313 Subscript colLow, Subscript colHi)
01314 {
01315 Subscript lowerBound = A_.lbound();
01316
01317 if ( rowHi < rowLow || rowLow < lowerBound || rowHi > A_.num_rows())
01318 throw RuntimeError("Bad row index range", "SCPPNT::const_Region2D<Array2D>::initialize()");
01319 if ( colHi < colLow || colLow < lowerBound || colHi > A_.num_columns())
01320 throw RuntimeError("Bad column index range", "SCPPNT::const_Region2D<Array2D>::initialize()");
01321
01322 offset_[0] = rowLow-lowerBound;
01323 offset_[1] = colLow-lowerBound;
01324 dim_[0] = rowHi-rowLow+1;
01325 dim_[1] = colHi-colLow+1;
01326
01327
01328 const_row_ = new const_row_iterator[dim_[0]+1];
01329
01330
01331 Subscript i;
01332 for (i=0; i<dim_[0]; i++)
01333 {
01334 const_row_[i] = A_.row(i+offset_[0]+1) + offset_[1];
01335 #ifdef SCPPNT_BOUNDS_CHECK
01336 const_row_[i].set_size(dim_[1]);
01337 #endif
01338 }
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352 const_column_ = new const_column_iterator[dim_[1]+1];
01353
01354
01355 for (i=0; i<dim_[1]; ++i)
01356 {
01357 const_column_[i] = A_.column(i+offset_[1]+1) + offset_[0];
01358 #ifdef SCPPNT_BOUNDS_CHECK
01359 const_column_[i].set_size(dim_[0]);
01360 #endif
01361 }
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373 }
01374
01375
01376 template <class Array2D>
01377 inline typename const_Region2D<Array2D>::const_rows_iterator const_Region2D<Array2D>::begin_rows() const
01378 {
01379 #ifdef SCPPNT_BOUNDS_CHECK
01380 return const_rows_iterator(const_row_, 1, dim_[0]);
01381 #else
01382 return const_row_;
01383 #endif
01384 }
01385
01386
01387 template <class Array2D>
01388 inline typename const_Region2D<Array2D>::const_columns_iterator const_Region2D<Array2D>::begin_columns() const
01389 {
01390 #ifdef SCPPNT_BOUNDS_CHECK
01391 return const_columns_iterator(const_column_, 1, dim_[1]);
01392 #else
01393 return const_column_;
01394 #endif
01395 }
01396
01397
01398 template <class Array2D>
01399 inline typename const_Region2D<Array2D>::const_rows_iterator const_Region2D<Array2D>::end_rows() const
01400 {
01401 #ifdef SCPPNT_BOUNDS_CHECK
01402 return const_rows_iterator(const_row_ + dim_[0], 1, 1);
01403 #else
01404 return (const_rows_iterator) row_ + dim_[0];
01405 #endif
01406 }
01407
01408
01409 template <class Array2D>
01410 inline typename const_Region2D<Array2D>::const_columns_iterator const_Region2D<Array2D>::end_columns() const
01411 {
01412 #ifdef SCPPNT_BOUNDS_CHECK
01413 return const_columns_iterator(const_column_ + dim_[1], 1, 1);
01414 #else
01415 return const_column_ + dim_[1];
01416 #endif
01417 }
01418
01419
01420 template <class Array2D>
01421 inline Subscript const_Region2D<Array2D>::diagonal_size(Subscript row, Subscript column) const
01422 {
01423 Subscript ncolumn = dim_[0] - column + 1;
01424 Subscript nrow = dim_[1] - row + 1;
01425
01426 return (ncolumn > nrow) ? nrow : ncolumn;
01427 }
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439 template <class Array2D>
01440 inline typename const_Region2D<Array2D>::const_diag_iterator
01441 const_Region2D<Array2D>::begin_diagonal(Subscript row, Subscript column) const
01442 {
01443 const_diag_iterator it = A_.begin_diagonal(row + offset_[0], column + offset_[1]);
01444
01445 #ifdef SCPPNT_BOUNDS_CHECK
01446 it.set_size(diagonal_size(row, column));
01447 #endif
01448
01449 return it;
01450 }
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463 template <class Array2D>
01464 inline typename const_Region2D<Array2D>::const_diag_iterator
01465 const_Region2D<Array2D>::end_diagonal(Subscript row, Subscript column) const
01466 {
01467
01468 const_diag_iterator it = begin_diagonal(row, column);
01469 return it + diagonal_size(row, column);
01470 }
01471
01472 #ifndef SCPPNT_NO_IO
01473
01474
01475
01476
01477
01478
01479
01480 template <class Array2D>
01481 std::ostream& operator<<(std::ostream &s, const const_Region2D<Array2D> &A)
01482 {
01483 Subscript M=A.num_rows();
01484 Subscript N=A.num_columns();
01485
01486 Subscript start = A.lbound();
01487 Subscript Mend = M + A.lbound() - 1;
01488 Subscript Nend = N + A.lbound() - 1;
01489
01490 s << M << " " << N << std::endl;
01491 for (Subscript i=start; i<=Mend; i++)
01492 {
01493 for (Subscript j=start; j<=Nend; j++)
01494 {
01495 s << A(i,j) << " ";
01496 }
01497 s << std::endl;
01498 }
01499
01500 return s;
01501 }
01502 #endif
01503
01504
01505
01506
01507 template <class Array2D, class M>
01508 inline M operator+(const const_Region2D<Array2D> &A,
01509 const M &B)
01510 {
01511 M sum(A.num_rows(), A.num_columns());
01512 matadd(A, B, sum);
01513 return sum;
01514 }
01515
01516
01517 template <class Array2D, class M>
01518 inline M operator-(const const_Region2D<Array2D> &A,
01519 const M &B)
01520 {
01521 M diff(A.num_rows(), A.num_columns());
01522 matsub(A, B, diff);
01523 return diff;
01524 }
01525
01526
01527 template <class Array2D, class M>
01528 inline M operator*(const const_Region2D<Array2D> &A,
01529 const M &B)
01530 {
01531 M c(A.num_rows(), B.num_columns());
01532 matmult(c, A, B);
01533 return c;
01534 }
01535
01536
01537 template <class T, class Array2D>
01538 inline Vector<T> operator*(const const_Region2D<Array2D> &A, const Vector<T> &x)
01539 {
01540 return matrix_times_vector< const_Region2D<Array2D>, Vector<T> >(A, x);
01541 }
01542
01543 }
01544
01545 #endif
01546