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
00044
00045
00046
00047 #ifndef SCPPNT_CMAT_H
00048 #define SCPPNT_CMAT_H
00049
00050 #ifdef SCPPNT_NO_DIR_PREFIX
00051 #include "scppnt.h"
00052 #include "slice_pointer.h"
00053 #include "scppnt_error.h"
00054 #include "matop.h"
00055 #include "vec.h"
00056 #else
00057 #include "scppnt/scppnt.h"
00058 #include "scppnt/slice_pointer.h"
00059 #include "scppnt/scppnt_error.h"
00060 #include "scppnt/matop.h"
00061 #include "scppnt/vec.h"
00062 #endif
00063
00064
00065
00066 #ifndef SCPPNT_NO_IO
00067 #include <string>
00068 #include <iostream>
00069 #include <strstream>
00070 #endif
00071
00072
00073
00074
00075
00076 #ifdef SCPPNT_USE_REGIONS
00077 #ifdef SCPPNT_NO_DIR_PREFIX
00078 #include "region2d.h"
00079 #else
00080 #include "scppnt/region2d.h"
00081 #endif
00082 #endif
00083
00084 namespace SCPPNT
00085 {
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 template<class T> class Matrix
00111 {
00112
00113 public:
00114
00115 typedef Subscript size_type;
00116 typedef T value_type;
00117 typedef T element_type;
00118 typedef T* pointer;
00119 typedef T& reference;
00120 typedef const T& const_reference;
00121
00122 #ifdef SCPPNT_BOUNDS_CHECK
00123 typedef slice_pointer_base<T, T*, T&> iterator;
00124 typedef slice_pointer_base<T, const T*, const T&> const_iterator;
00125 #else
00126 typedef T* iterator;
00127 typedef const T* const_iterator;
00128 #endif
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142 #ifdef SCPPNT_BOUNDS_CHECK
00143 typedef slice_pointer_base<T, T*, T&> row_iterator;
00144 typedef slice_pointer_base<T, const T*, const T&> const_row_iterator;
00145 #else
00146
00147 typedef T* row_iterator;
00148
00149
00150 typedef const T* const_row_iterator;
00151 #endif
00152
00153 typedef slice_pointer_base<T, T*, T&> column_iterator;
00154
00155
00156 typedef slice_pointer_base<T, const T*, const T&> const_column_iterator;
00157
00158
00159
00160
00161
00162 #ifdef SCPPNT_BOUNDS_CHECK
00163 typedef slice_pointer_base<row_iterator, row_iterator*, row_iterator&> rows_iterator;
00164 typedef slice_pointer_base<column_iterator, column_iterator*, column_iterator&>
00165 columns_iterator;
00166 typedef slice_pointer_base<const_row_iterator, const_row_iterator*, const_row_iterator&>
00167 const_rows_iterator;
00168 typedef slice_pointer_base<const_column_iterator, const_column_iterator*, const_column_iterator&>
00169 const_columns_iterator;
00170 #else
00171
00172 typedef row_iterator* rows_iterator;
00173
00174
00175 typedef const_row_iterator* const_rows_iterator;
00176
00177
00178 typedef column_iterator* columns_iterator;
00179
00180
00181 typedef const_column_iterator* const_columns_iterator;
00182 #endif
00183
00184
00185 typedef slice_pointer_base<T, T*, T&> diag_iterator;
00186
00187
00188 typedef slice_pointer_base<T, const T*, const T&> const_diag_iterator;
00189
00190
00191 Subscript lbound() const;
00192
00193
00194 Subscript size() const;
00195
00196
00197 Subscript dim(Subscript d) const;
00198
00199
00200 Subscript num_rows() const;
00201
00202
00203 Subscript num_columns() const;
00204
00205
00206
00207
00208 Matrix();
00209
00210
00211 Matrix(const Matrix<T> &A);
00212
00213
00214 Matrix(Subscript M, Subscript N);
00215
00216
00217 Matrix(Subscript M, Subscript N, const T value);
00218
00219 #ifndef SCPPNT_NO_IO
00220
00221 Matrix(Subscript M, Subscript N, const std::string &s);
00222 #endif
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246 template<class IT> Matrix(IT begin, Subscript number_rows, Subscript number_columns) :
00247 v_(0)
00248 {
00249
00250 initialize(number_rows, number_columns);
00251
00252 iterator i = this->begin();
00253 Subscript n = number_rows * number_columns;
00254 while (n--)
00255 {
00256 *i = *begin;
00257 ++begin;
00258 ++i;
00259 }
00260 }
00261
00262
00263 ~Matrix();
00264
00265
00266 Matrix<T>& newsize(Subscript M, Subscript N);
00267
00268
00269
00270
00271
00272 Matrix<T>& operator=(const Matrix<T> &A);
00273
00274
00275 Matrix<T>& operator=(const T& scalar);
00276
00277
00278
00279
00280
00281 row_iterator operator[](Subscript i);
00282
00283
00284 const_row_iterator operator[](Subscript i) const;
00285
00286
00287 reference operator()(Subscript i);
00288
00289
00290 const_reference operator()(Subscript i) const;
00291
00292
00293 reference operator()(Subscript i, Subscript j);
00294
00295
00296 const_reference operator()(Subscript i, Subscript j) const;
00297
00298
00299
00300
00301 rows_iterator begin_rows();
00302
00303
00304 const_rows_iterator begin_rows() const;
00305
00306
00307 columns_iterator begin_columns();
00308
00309
00310 const_columns_iterator begin_columns() const;
00311
00312
00313 rows_iterator end_rows();
00314
00315
00316 const_rows_iterator end_rows() const;
00317
00318
00319 columns_iterator end_columns();
00320
00321
00322 const_columns_iterator end_columns() const;
00323
00324
00325
00326
00327
00328
00329
00330
00331 row_iterator begin_row(Subscript index);
00332
00333
00334 const_row_iterator begin_row(Subscript index) const;
00335
00336
00337 column_iterator begin_column(Subscript index);
00338
00339
00340 const_column_iterator begin_column(Subscript index) const;
00341
00342
00343 row_iterator end_row(Subscript index);
00344
00345
00346 const_row_iterator end_row(Subscript index) const;
00347
00348
00349 column_iterator end_column(Subscript index);
00350
00351
00352 const_column_iterator end_column(Subscript index) const;
00353
00354
00355 iterator begin();
00356
00357
00358 const_iterator begin() const;
00359
00360
00361 iterator end();
00362
00363
00364 const_iterator end() const;
00365
00366
00367 diag_iterator begin_diagonal(Subscript row, Subscript column);
00368
00369
00370 const_diag_iterator begin_diagonal(Subscript row, Subscript column) const;
00371
00372
00373 diag_iterator end_diagonal(Subscript row, Subscript column);
00374
00375
00376 const_diag_iterator end_diagonal(Subscript row, Subscript column) const;
00377
00378 #ifndef SCPPNT_BOUNDS_CHECK
00379
00380
00381
00382
00383 operator T**()
00384 {
00385 return row_;
00386 }
00387
00388
00389 operator const T**() const
00390 {
00391 return (const T**) row_;
00392 }
00393 #endif
00394
00395
00396
00397
00398
00399
00400
00401
00402 template<class MAT> Matrix<T>& operator+=(const MAT &rhs)
00403 {
00404 matadd(*this, rhs);
00405 return *this;
00406 }
00407
00408
00409 template<class MAT> Matrix<T>& operator-=(const MAT &rhs)
00410 {
00411 matsub(*this, rhs);
00412 return *this;
00413 }
00414
00415
00416 template<class MAT> Matrix<T>& operator*=(const MAT &rhs)
00417 {
00418 matmult_assign(*this, rhs);
00419 return *this;
00420 }
00421
00422
00423 Matrix<T>& operator+=(const T& value);
00424
00425
00426 Matrix<T>& operator-=(const T& value);
00427
00428
00429 Matrix<T>& operator*=(const T& value);
00430
00431
00432 Matrix<T>& operator/=(const T& value);
00433
00434
00435
00436 #ifdef SCPPNT_USE_REGIONS
00437
00438
00439 typedef Region2D<Matrix<T> > Region;
00440
00441
00442 Region operator()(const Index1D &I, const Index1D &J)
00443 {
00444 return Region(*this, I,J);
00445 }
00446
00447
00448 typedef const_Region2D< Matrix<T> > const_Region;
00449
00450
00451 const_Region operator()(const Index1D &I, const Index1D &J) const
00452 {
00453 return const_Region(*this, I,J);
00454 }
00455
00456 #endif
00457
00458 private:
00459
00460 Subscript m_;
00461 Subscript n_;
00462 Subscript mn_;
00463 T* v_;
00464 T* vm1_;
00465 row_iterator *row_;
00466 T** rowm1_;
00467 column_iterator *column_;
00468 const_column_iterator *const_column_;
00469
00470 #ifdef SCPPNT_BOUNDS_CHECK
00471 const_row_iterator *const_row_;
00472 #endif
00473
00474
00475
00476
00477 void initialize(Subscript M, Subscript N);
00478
00479
00480 void copy(const T* v);
00481
00482
00483 void set(const T& val);
00484
00485
00486 void destroy();
00487
00488
00489 Subscript diagonal_size(Subscript row, Subscript column) const;
00490
00491 };
00492
00493 template<class T> inline Subscript Matrix<T>::lbound() const
00494 {
00495 return 1;
00496 }
00497
00498
00499
00500 template<class T> Matrix<T>& Matrix<T>::operator=(const Matrix<T> &A)
00501 {
00502 if (v_ == A.v_)
00503 return *this;
00504
00505 if (m_ == A.m_ && n_ == A.n_)
00506 copy(A.v_);
00507 else
00508 {
00509 destroy();
00510 initialize(A.m_, A.n_);
00511 copy(A.v_);
00512 }
00513
00514 return *this;
00515 }
00516
00517 template <class T>
00518 Matrix<T>& Matrix<T>::operator=(const T& scalar)
00519 {
00520 set(scalar);
00521 return *this;
00522 }
00523
00524 template <class T>
00525 inline Subscript Matrix<T>::size() const
00526 {
00527 return mn_;
00528 }
00529
00530 template <class T>
00531 Subscript Matrix<T>::dim(Subscript d) const
00532 {
00533 if (d < 1 || d > 2) throw InvalidArgument("argument must be 1 or 2", "SCPPNT::Matrix::dim");
00534 return (d==1) ? m_ : n_;
00535 }
00536
00537 template <class T>
00538 inline Subscript Matrix<T>::num_rows() const
00539 { return m_;}
00540 template <class T>
00541 inline Subscript Matrix<T>::num_columns() const
00542 { return n_;}
00543
00544 template <class T>
00545 void Matrix<T>::initialize(Subscript M, Subscript N)
00546 {
00547
00548 if (v_ != 0) throw LogicError(0, "v_ is not NULL", "SCPPNT::Matrix::initialize");
00549
00550 mn_ = M*N;
00551 m_ = M;
00552 n_ = N;
00553
00554
00555 v_ = new T[mn_];
00556 vm1_ = v_ - 1;
00557
00558
00559
00560 rowm1_ = new T*[M+1];
00561 row_ = new row_iterator[M+1];
00562 #ifdef SCPPNT_BOUNDS_CHECK
00563 const_row_ = new const_row_iterator[M+1];
00564 #endif
00565
00566
00567 T* p = v_;
00568 Subscript i;
00569 for (i=0; i<=M; i++)
00570 {
00571 #ifdef SCPPNT_BOUNDS_CHECK
00572 row_[i].set(p, 1, N);
00573 const_row_[i].set(p, 1, N);
00574 #else
00575 row_[i] = p;
00576 #endif
00577
00578 rowm1_[i] = p-1;
00579 p += N;
00580 }
00581
00582 rowm1_ --;
00583
00584
00585
00586 column_ = new column_iterator[N+1];
00587 const_column_ = new const_column_iterator[N+1];
00588
00589
00590 p = v_;
00591 for (i=0; i<=N; ++i)
00592 {
00593 column_[i].set(p, N, M);
00594 const_column_[i].set(p, N, M);
00595 ++p;
00596 }
00597
00598 }
00599
00600 template <class T>
00601 void Matrix<T>::destroy()
00602 {
00603
00604 if (v_ == 0) return;
00605
00606
00607 delete [] (v_);
00608 v_ = 0;
00609
00610 if (row_ != 0) delete [] (row_);
00611 if (column_ != 0) delete [] (column_);
00612
00613 #ifdef SCPPNT_BOUNDS_CHECK
00614 if (const_row_ != 0) delete [] (const_row_);
00615 #endif
00616 if (const_column_ != 0) delete [] (const_column_);
00617
00618 if (rowm1_ != 0)
00619 {
00620 rowm1_ ++;
00621 if (rowm1_ != 0 ) delete [] (rowm1_);
00622 }
00623
00624 }
00625
00626 template <class T>
00627 Matrix<T>& Matrix<T>::newsize(Subscript M, Subscript N)
00628 {
00629 if (num_rows() == M && num_columns() == N)
00630 return *this;
00631
00632 destroy();
00633 initialize(M,N);
00634
00635 return *this;
00636 }
00637
00638 template <class T>
00639 void Matrix<T>::copy(const T* v)
00640 {
00641 Subscript N = m_ * n_;
00642 Subscript i;
00643
00644 #ifdef SCPPNT_UNROLL_LOOPS
00645 Subscript Nmod4 = N & 3;
00646 Subscript N4 = N - Nmod4;
00647
00648 for (i=0; i<N4; i+=4)
00649 {
00650 v_[i] = v[i];
00651 v_[i+1] = v[i+1];
00652 v_[i+2] = v[i+2];
00653 v_[i+3] = v[i+3];
00654 }
00655
00656 for (i=N4; i< N; i++)
00657 v_[i] = v[i];
00658 #else
00659
00660 for (i=0; i< N; i++)
00661 v_[i] = v[i];
00662 #endif
00663 }
00664
00665 template <class T>
00666 void Matrix<T>::set(const T& val)
00667 {
00668 Subscript N = m_ * n_;
00669 Subscript i;
00670
00671 #ifdef SCPPNT_UNROLL_LOOPS
00672 Subscript Nmod4 = N & 3;
00673 Subscript N4 = N - Nmod4;
00674
00675 for (i=0; i<N4; i+=4)
00676 {
00677 v_[i] = val;
00678 v_[i+1] = val;
00679 v_[i+2] = val;
00680 v_[i+3] = val;
00681 }
00682
00683 for (i=N4; i< N; i++)
00684 v_[i] = val;
00685 #else
00686 iterator pi = begin();
00687 for (i=N; i--; ++pi)
00688 *pi = val;
00689
00690 #endif
00691 }
00692
00693
00694
00695 template <class T>
00696 Matrix<T>::Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0), column_(0), const_column_(0)
00697 {
00698 #ifdef SCPPNT_BOUNDS_CHECK
00699 const_row_ = 0;
00700 #endif
00701 }
00702
00703 template <class T>
00704 Matrix<T>::Matrix(const Matrix<T> &A) : v_(0)
00705 {
00706 initialize(A.m_, A.n_);
00707 copy(A.v_);
00708 }
00709
00710 template <class T>
00711 Matrix<T>::Matrix(Subscript M, Subscript N) : v_(0)
00712 {
00713 initialize(M,N);
00714 }
00715
00716 template <class T>
00717 Matrix<T>::Matrix(Subscript M, Subscript N, const T value) : v_(0)
00718 {
00719 initialize(M,N);
00720 set(value);
00721 }
00722
00723 #ifndef SCPPNT_NO_IO
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735 template <class T>
00736 Matrix<T>::Matrix(Subscript M, Subscript N, const std::string &s) : v_(0)
00737 {
00738 initialize(M,N);
00739 std::istrstream ins(s.c_str());
00740
00741 Subscript i, j;
00742
00743 for (i=0; i<M; i++)
00744 for (j=0; j<N; j++)
00745 ins >> row_[i][j];
00746 }
00747 #endif
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759 template <class T>
00760 Matrix<T>::~Matrix()
00761 {
00762 destroy();
00763 }
00764
00765
00766 template <class T>
00767 inline typename Matrix<T>::row_iterator Matrix<T>::operator[](Subscript i)
00768 {
00769 #ifdef SCPPNT_BOUNDS_CHECK
00770 if (i < 0 || i >= m_) throw BoundsError("SCPPNT::Matrix::operator[]");
00771 #endif
00772 return row_[i];
00773 }
00774
00775 template <class T>
00776 inline typename Matrix<T>::const_row_iterator Matrix<T>::operator[](Subscript i) const
00777 {
00778 #ifdef SCPPNT_BOUNDS_CHECK
00779 if (i < 0 || i >= m_) throw BoundsError("SCPPNT::Matrix::operator[] const");
00780 #endif
00781 return row_[i];
00782 }
00783
00784 template <class T>
00785 inline typename Matrix<T>::reference Matrix<T>::operator()(Subscript i)
00786 {
00787 #ifdef SCPPNT_BOUNDS_CHECK
00788 if (i < 1 || i > mn_) throw BoundsError("SCPPNT::Matrix::operator(i)");
00789 #endif
00790 return vm1_[i];
00791 }
00792
00793 template <class T>
00794 inline typename Matrix<T>::const_reference Matrix<T>::operator()(Subscript i) const
00795 {
00796 #ifdef SCPPNT_BOUNDS_CHECK
00797 if (i < 1 || i > mn_) throw BoundsError("SCPPNT::Matrix::operator(i) const");
00798 #endif
00799 return vm1_[i];
00800 }
00801
00802 template <class T>
00803 inline typename Matrix<T>::reference Matrix<T>::operator()(Subscript i, Subscript j)
00804 {
00805 #ifdef SCPPNT_BOUNDS_CHECK
00806 if (i < 1 || i > m_ || j < 1 || j > n_) throw BoundsError("SCPPNT::Matrix::operator(i,j)");
00807 #endif
00808 return rowm1_[i][j];
00809 }
00810
00811 template <class T>
00812 inline typename Matrix<T>::const_reference Matrix<T>::operator() (Subscript i, Subscript j) const
00813 {
00814 #ifdef SCPPNT_BOUNDS_CHECK
00815 if (i < 1 || i > m_ || j < 1 || j > n_) throw BoundsError("SCPPNT::Matrix::operator(i,j) const");
00816 #endif
00817 return rowm1_[i][j];
00818 }
00819
00820
00821
00822
00823 template <class T>
00824 inline typename Matrix<T>::iterator Matrix<T>::begin()
00825 {
00826 #ifdef SCPPNT_BOUNDS_CHECK
00827 return iterator(v_, 1, mn_);
00828 #else
00829 return v_;
00830 #endif
00831 }
00832
00833 template <class T>
00834 inline typename Matrix<T>::const_iterator Matrix<T>::begin() const
00835 {
00836 #ifdef SCPPNT_BOUNDS_CHECK
00837 return const_iterator(v_, 1, mn_);
00838 #else
00839 return v_;
00840 #endif
00841 }
00842
00843 template <class T>
00844 inline typename Matrix<T>::iterator Matrix<T>::end()
00845 {
00846 #ifdef SCPPNT_BOUNDS_CHECK
00847
00848
00849
00850 return iterator(v_, 1, mn_) += mn_;
00851 #else
00852 return v_ + mn_;
00853 #endif
00854 }
00855
00856 template <class T>
00857 inline typename Matrix<T>::const_iterator Matrix<T>::end() const
00858 {
00859 #ifdef SCPPNT_BOUNDS_CHECK
00860
00861
00862
00863 return const_iterator(v_, 1, mn_) += mn_;
00864 #else
00865 return v_ + mn_;
00866 #endif
00867 }
00868
00869
00870 template <class T>
00871 inline typename Matrix<T>::rows_iterator Matrix<T>::begin_rows()
00872 {
00873 #ifdef SCPPNT_BOUNDS_CHECK
00874 return rows_iterator(row_, 1, m_);
00875 #else
00876 return row_;
00877 #endif
00878 }
00879
00880 template <class T>
00881 inline typename Matrix<T>::const_rows_iterator Matrix<T>::begin_rows() const
00882 {
00883 #ifdef SCPPNT_BOUNDS_CHECK
00884 return const_rows_iterator(const_row_, 1, m_);
00885 #else
00886 return (const_rows_iterator) row_;
00887 #endif
00888 }
00889
00890 template <class T>
00891 inline typename Matrix<T>::columns_iterator Matrix<T>::begin_columns()
00892 {
00893 #ifdef SCPPNT_BOUNDS_CHECK
00894 return columns_iterator(column_, 1, n_);
00895 #else
00896 return column_;
00897 #endif
00898 }
00899
00900 template <class T>
00901 inline typename Matrix<T>::const_columns_iterator Matrix<T>::begin_columns() const
00902 {
00903 #ifdef SCPPNT_BOUNDS_CHECK
00904 return const_columns_iterator(const_column_, 1, n_);
00905 #else
00906 return const_column_;
00907 #endif
00908 }
00909
00910 template <class T>
00911 inline typename Matrix<T>::rows_iterator Matrix<T>::end_rows()
00912 {
00913 #ifdef SCPPNT_BOUNDS_CHECK
00914 return rows_iterator(row_ + m_, 1, 1);
00915 #else
00916 return row_ + m_;
00917 #endif
00918 }
00919
00920 template <class T>
00921 inline typename Matrix<T>::const_rows_iterator Matrix<T>::end_rows() const
00922 {
00923 #ifdef SCPPNT_BOUNDS_CHECK
00924 return const_rows_iterator(const_row_ + m_, 1, 1);
00925 #else
00926 return (const_rows_iterator) row_ + m_;
00927 #endif
00928 }
00929
00930 template <class T>
00931 inline typename Matrix<T>::columns_iterator Matrix<T>::end_columns()
00932 {
00933 #ifdef SCPPNT_BOUNDS_CHECK
00934 return columns_iterator(column_ + n_, 1, 1);
00935 #else
00936 return column_ + n_;
00937 #endif
00938 }
00939
00940 template <class T>
00941 inline typename Matrix<T>::const_columns_iterator Matrix<T>::end_columns() const
00942 {
00943 #ifdef SCPPNT_BOUNDS_CHECK
00944 return const_columns_iterator(const_column_ + n_, 1, 1);
00945 #else
00946 return const_column_ + n_;
00947 #endif
00948 }
00949
00950 template <class T>
00951 inline typename Matrix<T>::row_iterator Matrix<T>::begin_row(Subscript index)
00952 {
00953 #ifdef SCPPNT_BOUNDS_CHECK
00954 if (index < 1 || index > m_) throw BoundsError("SCPPNT::Matrix::begin_row()");
00955 #endif
00956 return row_[index-1];
00957 }
00958
00959 template <class T>
00960 inline typename Matrix<T>::const_row_iterator Matrix<T>::begin_row(Subscript index) const
00961 {
00962 #ifdef SCPPNT_BOUNDS_CHECK
00963 if (index < 1 || index > m_) throw BoundsError("SCPPNT::Matrix::begin_row() const");
00964 return const_row_[index-1];
00965 #else
00966 return (const_row_iterator) row_[index-1];
00967 #endif
00968 }
00969
00970 template <class T>
00971 inline typename Matrix<T>::column_iterator Matrix<T>::begin_column(Subscript index)
00972 {
00973 #ifdef SCPPNT_BOUNDS_CHECK
00974 if (index < 1 || index > n_) throw BoundsError("SCPPNT::Matrix::begin_column()");
00975 #endif
00976 return column_[index-1];
00977 }
00978
00979 template <class T>
00980 inline typename Matrix<T>::const_column_iterator Matrix<T>::begin_column(Subscript index) const
00981 {
00982 #ifdef SCPPNT_BOUNDS_CHECK
00983 if (index < 1 || index > n_) throw BoundsError("SCPPNT::Matrix::begin_column() const");
00984 #endif
00985 return const_column_[index-1];
00986 }
00987
00988 template <class T>
00989 inline typename Matrix<T>::row_iterator Matrix<T>::end_row(Subscript index)
00990 {
00991 #ifdef SCPPNT_BOUNDS_CHECK
00992 if (index < 1 || index > m_) throw BoundsError("SCPPNT::Matrix::end_row()");
00993 #endif
00994 return row_[index];
00995 }
00996
00997 template <class T>
00998 inline typename Matrix<T>::const_row_iterator Matrix<T>::end_row(Subscript index) const
00999 {
01000 #ifdef SCPPNT_BOUNDS_CHECK
01001 if (index < 1 || index > m_) throw BoundsError("SCPPNT::Matrix::end_row() const");
01002 return const_row_[index];
01003 #else
01004 return (const_row_iterator) row_[index];
01005 #endif
01006 }
01007
01008 template <class T>
01009 inline typename Matrix<T>::column_iterator Matrix<T>::end_column(Subscript index)
01010 {
01011 #ifdef SCPPNT_BOUNDS_CHECK
01012 if (index < 1 || index > n_) throw BoundsError("SCPPNT::Matrix::end_column()");
01013 #endif
01014 return column_[index];
01015 }
01016
01017 template <class T>
01018 inline typename Matrix<T>::const_column_iterator Matrix<T>::end_column(Subscript index) const
01019 {
01020 #ifdef SCPPNT_BOUNDS_CHECK
01021 if (index < 1 || index > n_) throw BoundsError("SCPPNT::Matrix::end_column() const");
01022 #endif
01023 return const_column_[index];
01024 }
01025
01026
01027 template <class T>
01028 inline Subscript Matrix<T>::diagonal_size(Subscript row, Subscript column) const
01029 {
01030 Subscript ncolumn = n_ - column + 1;
01031 Subscript nrow = m_ - row + 1;
01032
01033 return (ncolumn > nrow) ? nrow : ncolumn;
01034 }
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046 template <class T>
01047 inline typename Matrix<T>::diag_iterator Matrix<T>::begin_diagonal(Subscript row, Subscript column)
01048 {
01049 Subscript num_elem = 0;
01050
01051 #ifdef SCPPNT_BOUNDS_CHECK
01052
01053 if (row < 1 || row > m_ || column < 1 || column > n_) BoundsError("SCPPNT::Matrix::begin_diagonal(row, column)");
01054 num_elem = diagonal_size(row, column);
01055
01056 #endif
01057
01058
01059
01060 return diag_iterator(v_ + column-1 + n_*(row-1), n_ + 1, num_elem);
01061 }
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073 template <class T>
01074 inline typename Matrix<T>::const_diag_iterator Matrix<T>::begin_diagonal(Subscript row, Subscript column) const
01075 {
01076 Subscript num_elem = 0;
01077
01078 #ifdef SCPPNT_BOUNDS_CHECK
01079
01080 if (row < 1 || row > m_ || column < 1 || column > n_) BoundsError("SCPPNT::Matrix::begin_diagonal(row, column)");
01081 num_elem = diagonal_size(row, column);
01082
01083 #endif
01084
01085
01086
01087 return const_diag_iterator(v_ + column-1 + n_*(row-1), n_ + 1, num_elem);
01088 }
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101 template <class T>
01102 inline typename Matrix<T>::diag_iterator Matrix<T>::end_diagonal(Subscript row, Subscript column)
01103 {
01104 diag_iterator it = begin_diagonal(row, column);
01105 return it + diagonal_size(row, column);
01106 }
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118 template <class T>
01119 inline typename Matrix<T>::const_diag_iterator Matrix<T>::end_diagonal(Subscript row, Subscript column) const
01120 {
01121 const_diag_iterator it = begin_diagonal(row, column);
01122 return it + diagonal_size(row, column);
01123 }
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146 template <class T>
01147 Matrix<T>& Matrix<T>::operator+=(const T& value)
01148 {
01149 double *thisiter = v_;
01150 for (int i = mn_; i--; ++thisiter)
01151 {
01152 *thisiter += value;
01153 }
01154 return *this;
01155 }
01156
01157
01158 template <class T>
01159 Matrix<T>& Matrix<T>::operator-=(const T& value)
01160 {
01161 double *thisiter = v_;
01162 for (int i = mn_; i--; ++thisiter)
01163 {
01164 *thisiter -= value;
01165 }
01166 return *this;
01167 }
01168
01169
01170 template <class T>
01171 Matrix<T>& Matrix<T>::operator*=(const T& value)
01172 {
01173 double *thisiter = v_;
01174 for (Subscript i = mn_; i--; ++thisiter)
01175 {
01176 *thisiter *= value;
01177 }
01178 return *this;
01179 }
01180
01181
01182 template <class T>
01183 Matrix<T>& Matrix<T>::operator/=(const T& value)
01184 {
01185 double *thisiter = v_;
01186 for (Subscript i = mn_; i--; ++thisiter)
01187 {
01188 *thisiter /= value;
01189 }
01190 return *this;
01191 }
01192
01193 #ifndef SCPPNT_NO_IO
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203 template <class T>
01204 std::ostream& operator<<(std::ostream &s, const Matrix<T> &A)
01205 {
01206 Subscript M = A.num_rows();
01207 Subscript N = A.num_columns();
01208
01209 s << M << " " << N << "\n";
01210
01211 for (Subscript i=0; i<M; i++)
01212 {
01213 for (Subscript j=0; j<N; j++)
01214 {
01215 s << A[i][j] << " ";
01216 }
01217 s << "\n";
01218 }
01219 return s;
01220 }
01221
01222
01223
01224
01225
01226
01227
01228 template <class T>
01229 std::istream& operator>>(std::istream &s, Matrix<T> &A)
01230 {
01231 Subscript M, N;
01232
01233 s >> M >> N;
01234
01235 A.newsize(M,N);
01236
01237 for (Subscript i=0; i<M; i++)
01238 for (Subscript j=0; j<N; j++)
01239 {
01240 s >> A[i][j];
01241 }
01242 return s;
01243 }
01244 #endif // SCPPNT_NO_IO
01245
01246
01247
01248 template <class T, class M>
01249 inline Matrix<T> operator+(const Matrix<T> &A, const M &B)
01250 {
01251 return Matrix<T>(A) += B;
01252 }
01253
01254
01255 template <class T, class M>
01256 inline Matrix<T> operator-(const Matrix<T> &A, const M &B)
01257 {
01258 return Matrix<T>(A) -= B;
01259 }
01260
01261
01262 template <class T, class M>
01263 inline Matrix<T> operator*(const Matrix<T> &A, const M &B)
01264 {
01265 return matmult<Matrix<T>, Matrix<T>, M>(A,B);
01266 }
01267
01268
01269 template <class T>
01270 inline Vector<T> operator*(const Matrix<T> &A, const Vector<T> &x)
01271 {
01272 return matrix_times_vector< Matrix<T>, Vector<T> >(A, x);
01273 }
01274
01275 }
01276
01277 #endif
01278