00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00094 #ifndef matpackI_h
00095 #define matpackI_h
00096
00097 #include <iostream>
00098 #include <iomanip>
00099 #include "arts.h"
00100
00111 class Joker {
00112
00113 };
00114
00115
00116 extern Joker joker;
00117
00118
00119 class VectorView;
00120
00121
00122 class SparseMatrixView;
00123
00139 class Range{
00140 public:
00141
00142 Range(Index start, Index extent, Index stride=1);
00143 Range(Index start, Joker joker, Index stride=1);
00144 Range(Joker joker, Index stride=1);
00145 Range(Index max_size, const Range& r);
00146 Range(const Range& p, const Range& n);
00147
00148
00149 friend class ConstVectorView;
00150 friend class VectorView;
00151 friend class Vector;
00152 friend class ConstMatrixView;
00153 friend class MatrixView;
00154 friend class Matrix;
00155 friend class Iterator2D;
00156 friend class Iterator3D;
00157 friend class ConstIterator2D;
00158 friend class ConstIterator3D;
00159 friend class SparseMatrixView;
00160 friend class ConstTensor3View;
00161 friend class Tensor3View;
00162 friend class Tensor3;
00163 friend std::ostream& operator<<(std::ostream& os, const SparseMatrixView& v);
00164
00165
00166 private:
00168 Index mstart;
00171 Index mextent;
00173 Index mstride;
00174 };
00175
00178 class Iterator1D {
00179 public:
00180
00181 Iterator1D();
00182 Iterator1D(const Iterator1D& o);
00183 Iterator1D(Numeric *x, Index stride);
00184
00185
00186 Iterator1D& operator++();
00187 Numeric& operator*() const;
00188 bool operator!=(const Iterator1D& other) const;
00189
00190 private:
00192 Numeric *mx;
00194 Index mstride;
00195 };
00196
00199 class ConstIterator1D {
00200 public:
00201
00202 ConstIterator1D();
00203 ConstIterator1D(const ConstIterator1D& o);
00204 ConstIterator1D(Numeric *x, Index stride);
00205
00206
00207 ConstIterator1D& operator++();
00208 const Numeric& operator*() const;
00209 bool operator!=(const ConstIterator1D& other) const;
00210
00211 private:
00213 const Numeric *mx;
00215 Index mstride;
00216 };
00217
00218
00219 class Vector;
00220
00221
00222 class MatrixView;
00223
00224
00225 class ConstMatrixView;
00226
00232 class ConstVectorView {
00233 public:
00234
00235
00236 Index nelem() const;
00237 Numeric sum() const;
00238
00239
00240 Numeric operator[](Index n) const;
00241 ConstVectorView operator[](const Range& r) const;
00242
00243
00244 ConstIterator1D begin() const;
00245 ConstIterator1D end() const;
00246
00247
00248 operator ConstMatrixView() const;
00249
00250
00251 friend class VectorView;
00252 friend class ConstIterator2D;
00253 friend class ConstMatrixView;
00254 friend class ConstTensor3View;
00255
00256 protected:
00257
00258 ConstVectorView();
00259 ConstVectorView(Numeric *data, const Range& range);
00260 ConstVectorView(Numeric *data, const Range& p, const Range& n);
00261
00262
00263
00265 Range mrange;
00267 Numeric *mdata;
00268 };
00269
00281 class VectorView : public ConstVectorView {
00282 public:
00283
00284
00285 Numeric operator[](Index n) const;
00286 ConstVectorView operator[](const Range& r) const;
00287
00288 Numeric& operator[](Index n);
00289 VectorView operator[](const Range& r);
00290
00291
00292 ConstIterator1D begin() const;
00293 ConstIterator1D end() const;
00294
00295 Iterator1D begin();
00296 Iterator1D end();
00297
00298
00299 VectorView operator=(const ConstVectorView& v);
00300 VectorView operator=(const VectorView& v);
00301 VectorView operator=(const Vector& v);
00302 VectorView operator=(const Array<Numeric>& v);
00303 VectorView operator=(Numeric x);
00304
00305
00306 VectorView operator*=(Numeric x);
00307 VectorView operator/=(Numeric x);
00308 VectorView operator+=(Numeric x);
00309 VectorView operator-=(Numeric x);
00310
00311 VectorView operator*=(const ConstVectorView& x);
00312 VectorView operator/=(const ConstVectorView& x);
00313 VectorView operator+=(const ConstVectorView& x);
00314 VectorView operator-=(const ConstVectorView& x);
00315
00316
00317 operator MatrixView();
00318
00319
00320 friend class ConstIterator2D;
00321 friend class Iterator2D;
00322 friend class MatrixView;
00323 friend class Tensor3View;
00324
00325 protected:
00326
00327 VectorView();
00328 VectorView(Numeric *data, const Range& range);
00329 VectorView(Numeric *data, const Range& p, const Range& n);
00330 };
00331
00335 class Iterator2D {
00336 public:
00337
00338 Iterator2D();
00339 Iterator2D(const Iterator2D& o);
00340 Iterator2D(const VectorView& x, Index stride);
00341
00342
00343 Iterator2D& operator++();
00344 bool operator!=(const Iterator2D& other) const;
00345 VectorView* const operator->();
00346 VectorView& operator*();
00347
00348 private:
00350 VectorView msv;
00352 Index mstride;
00353 };
00354
00358 class ConstIterator2D {
00359 public:
00360
00361 ConstIterator2D();
00362 ConstIterator2D(const ConstIterator2D& o);
00363 ConstIterator2D(const ConstVectorView& x, Index stride);
00364
00365
00366 ConstIterator2D& operator++();
00367 bool operator!=(const ConstIterator2D& other) const;
00368 const ConstVectorView* operator->() const;
00369 const ConstVectorView& operator*() const;
00370
00371 private:
00373 ConstVectorView msv;
00375 Index mstride;
00376 };
00377
00389 class Vector : public VectorView {
00390 public:
00391
00392 Vector();
00393 explicit Vector(Index n);
00394 Vector(Index n, Numeric fill);
00395 Vector(Numeric start, Index extent, Numeric stride);
00396 Vector(const ConstVectorView& v);
00397 Vector(const Vector& v);
00398
00399
00400
00401 Vector& operator=(const Vector& v);
00402 Vector& operator=(const Array<Numeric>& v);
00403 Vector& operator=(Numeric x);
00404
00405
00406 void resize(Index n);
00407
00408
00409 ~Vector();
00410 };
00411
00412
00413 class Matrix;
00414
00415
00426 class ConstMatrixView {
00427 public:
00428
00429 Index nrows() const;
00430 Index ncols() const;
00431
00432
00433 Numeric operator()(Index r, Index c) const;
00434 ConstMatrixView operator()(const Range& r, const Range& c) const;
00435 ConstVectorView operator()(const Range& r, Index c) const;
00436 ConstVectorView operator()(Index r, const Range& c) const;
00437
00438
00439 ConstIterator2D begin() const;
00440 ConstIterator2D end() const;
00441
00442
00443 friend class ConstVectorView;
00444 friend class MatrixView;
00445 friend class ConstIterator3D;
00446 friend class ConstTensor3View;
00447 friend ConstMatrixView transpose(ConstMatrixView m);
00448
00449 protected:
00450
00451 ConstMatrixView();
00452 ConstMatrixView(Numeric *data, const Range& r, const Range& c);
00453 ConstMatrixView(Numeric *data,
00454 const Range& pr, const Range& pc,
00455 const Range& nr, const Range& nc);
00456
00457
00458
00460 Range mrr;
00462 Range mcr;
00464 Numeric *mdata;
00465 };
00466
00476 class MatrixView : public ConstMatrixView {
00477 public:
00478
00479
00480 Numeric operator()(Index r, Index c) const;
00481 ConstMatrixView operator()(const Range& r, const Range& c) const;
00482 ConstVectorView operator()(const Range& r, Index c) const;
00483 ConstVectorView operator()(Index r, const Range& c) const;
00484
00485 Numeric& operator()(Index r, Index c);
00486 MatrixView operator()(const Range& r, const Range& c);
00487 VectorView operator()(const Range& r, Index c);
00488 VectorView operator()(Index r, const Range& c);
00489
00490
00491 ConstIterator2D begin() const;
00492 ConstIterator2D end() const;
00493
00494 Iterator2D begin();
00495 Iterator2D end();
00496
00497
00498 MatrixView& operator=(const ConstMatrixView& v);
00499 MatrixView& operator=(const MatrixView& v);
00500 MatrixView& operator=(const Matrix& v);
00501 MatrixView& operator=(const ConstVectorView& v);
00502 MatrixView& operator=(Numeric x);
00503
00504
00505 MatrixView& operator*=(Numeric x);
00506 MatrixView& operator/=(Numeric x);
00507 MatrixView& operator+=(Numeric x);
00508 MatrixView& operator-=(Numeric x);
00509
00510 MatrixView& operator*=(const ConstMatrixView& x);
00511 MatrixView& operator/=(const ConstMatrixView& x);
00512 MatrixView& operator+=(const ConstMatrixView& x);
00513 MatrixView& operator-=(const ConstMatrixView& x);
00514
00515 MatrixView& operator*=(const ConstVectorView& x);
00516 MatrixView& operator/=(const ConstVectorView& x);
00517 MatrixView& operator+=(const ConstVectorView& x);
00518 MatrixView& operator-=(const ConstVectorView& x);
00519
00520
00521 friend class VectorView;
00522 friend class Iterator3D;
00523 friend class Tensor3View;
00524 friend ConstMatrixView transpose(ConstMatrixView m);
00525 friend MatrixView transpose(MatrixView m);
00526
00527 protected:
00528
00529 MatrixView();
00530 MatrixView(Numeric *data, const Range& r, const Range& c);
00531 MatrixView(Numeric *data,
00532 const Range& pr, const Range& pc,
00533 const Range& nr, const Range& nc);
00534 };
00535
00544 class Matrix : public MatrixView {
00545 public:
00546
00547 Matrix();
00548 Matrix(Index r, Index c);
00549 Matrix(Index r, Index c, Numeric fill);
00550 Matrix(const ConstMatrixView& v);
00551 Matrix(const Matrix& v);
00552
00553
00554 Matrix& operator=(const Matrix& x);
00555 Matrix& operator=(Numeric x);
00556 Matrix& operator=(const ConstVectorView& v);
00557
00558
00559 void resize(Index r, Index c);
00560
00561
00562 ~Matrix();
00563 };
00564
00565
00566
00567
00568
00569 inline void copy(ConstIterator1D origin,
00570 const ConstIterator1D& end,
00571 Iterator1D target);
00572
00573 inline void copy(Numeric x,
00574 Iterator1D target,
00575 const Iterator1D& end);
00576
00577 inline void copy(ConstIterator2D origin,
00578 const ConstIterator2D& end,
00579 Iterator2D target);
00580
00581 inline void copy(Numeric x,
00582 Iterator2D target,
00583 const Iterator2D& end);
00584
00585
00586
00587
00588 template<class base>
00589 class Array;
00590
00592 typedef Array<Vector> ArrayOfVector;
00593
00595 typedef Array<Matrix> ArrayOfMatrix;
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611 inline Range::Range(Index start, Index extent, Index stride) :
00612 mstart(start), mextent(extent), mstride(stride)
00613 {
00614
00615 assert( 0<=mstart );
00616
00617
00618
00619 assert( 0<=mextent );
00620
00621
00622
00623 }
00624
00627 inline Range::Range(Index start, Joker , Index stride) :
00628 mstart(start), mextent(-1), mstride(stride)
00629 {
00630
00631 assert( 0<=mstart );
00632 }
00633
00637 inline Range::Range(Joker , Index stride) :
00638 mstart(0), mextent(-1), mstride(stride)
00639 {
00640
00641 }
00642
00648 inline Range::Range(Index max_size, const Range& r) :
00649 mstart(r.mstart),
00650 mextent(r.mextent),
00651 mstride(r.mstride)
00652 {
00653
00654 assert( 0<=mstart );
00655
00656 assert( mstart<max_size );
00657
00658
00659 assert( 0!=mstride);
00660
00661
00662 if ( mextent<0 )
00663 {
00664 if ( 0<mstride )
00665 mextent = 1 + (max_size-1-mstart)/mstride;
00666 else
00667 mextent = 1 + (0-mstart)/mstride;
00668 }
00669 else
00670 {
00671 #ifndef NDEBUG
00672
00673 Index fin = mstart+(mextent-1)*mstride;
00674 assert( 0 <= fin );
00675 assert( fin < max_size );
00676 #endif
00677 }
00678 }
00679
00685 inline Range::Range(const Range& p, const Range& n) :
00686 mstart(p.mstart + n.mstart*p.mstride),
00687 mextent(n.mextent),
00688 mstride(p.mstride*n.mstride)
00689 {
00690
00691
00692
00693
00694
00695
00696
00697 Index prev_fin = p.mstart+(p.mextent-1)*p.mstride;
00698
00699
00700 assert( p.mstart<=mstart );
00701
00702 assert( mstart<=prev_fin );
00703
00704
00705 assert( 0!=mstride);
00706
00707
00708 if ( mextent<0 )
00709 {
00710 if ( 0<mstride )
00711 mextent = 1 + (prev_fin-mstart)/mstride;
00712 else
00713 mextent = 1 + (p.mstart-mstart)/mstride;
00714 }
00715 else
00716 {
00717 #ifndef NDEBUG
00718
00719 Index fin = mstart+(mextent-1)*mstride;
00720 assert( p.mstart <= fin );
00721 assert( fin <= prev_fin );
00722 #endif
00723 }
00724
00725 }
00726
00727
00728
00729
00730
00732 inline Iterator1D::Iterator1D()
00733 {
00734
00735 }
00736
00738 inline Iterator1D::Iterator1D(const Iterator1D& o) :
00739 mx(o.mx), mstride(o.mstride)
00740 {
00741
00742 }
00743
00745 inline Iterator1D::Iterator1D(Numeric *x, Index stride) :
00746 mx(x), mstride(stride)
00747 {
00748
00749 }
00750
00752 inline Iterator1D& Iterator1D::operator++()
00753 {
00754 mx += mstride;
00755 return *this;
00756 }
00757
00759 inline bool Iterator1D::operator!=(const Iterator1D& other) const
00760 {
00761 if (mx!=other.mx)
00762 return true;
00763 else
00764 return false;
00765 }
00766
00768 inline Numeric& Iterator1D::operator*() const
00769 {
00770 return *mx;
00771 }
00772
00773
00774
00775
00777 inline ConstIterator1D::ConstIterator1D()
00778 {
00779
00780 }
00781
00783 inline ConstIterator1D::ConstIterator1D(const ConstIterator1D& o) :
00784 mx(o.mx), mstride(o.mstride)
00785 {
00786
00787 }
00788
00790 inline ConstIterator1D::ConstIterator1D(Numeric *x, Index stride) :
00791 mx(x), mstride(stride)
00792 {
00793
00794 }
00795
00797 inline ConstIterator1D& ConstIterator1D::operator++()
00798 {
00799 mx += mstride;
00800 return *this;
00801 }
00802
00804 inline bool ConstIterator1D::operator!=(const ConstIterator1D& other) const
00805 {
00806 if (mx!=other.mx)
00807 return true;
00808 else
00809 return false;
00810 }
00811
00813 inline const Numeric& ConstIterator1D::operator*() const
00814 {
00815 return *mx;
00816 }
00817
00818
00819
00820
00822 inline Iterator2D::Iterator2D()
00823 {
00824
00825 }
00826
00828 inline Iterator2D::Iterator2D(const Iterator2D& o) :
00829 msv(o.msv), mstride(o.mstride)
00830 {
00831
00832 }
00833
00835 inline Iterator2D::Iterator2D(const VectorView& x, Index stride) :
00836 msv(x), mstride(stride)
00837 {
00838
00839 }
00840
00842 inline Iterator2D& Iterator2D::operator++()
00843 {
00844 msv.mdata += mstride;
00845 return *this;
00846 }
00847
00849 inline bool Iterator2D::operator!=(const Iterator2D& other) const
00850 {
00851 if ( msv.mdata + msv.mrange.mstart !=
00852 other.msv.mdata + other.msv.mrange.mstart )
00853 return true;
00854 else
00855 return false;
00856 }
00857
00860 inline VectorView* const Iterator2D::operator->()
00861 {
00862 return &msv;
00863 }
00864
00866 inline VectorView& Iterator2D::operator*()
00867 {
00868 return msv;
00869 }
00870
00871
00872
00873
00875 inline ConstIterator2D::ConstIterator2D()
00876 {
00877
00878 }
00879
00881 inline ConstIterator2D::ConstIterator2D(const ConstIterator2D& o) :
00882 msv(o.msv), mstride(o.mstride)
00883 {
00884
00885 }
00886
00888 inline ConstIterator2D::ConstIterator2D(const ConstVectorView& x, Index stride) :
00889 msv(x), mstride(stride)
00890 {
00891
00892 }
00893
00895 inline ConstIterator2D& ConstIterator2D::operator++()
00896 {
00897 msv.mdata += mstride;
00898 return *this;
00899 }
00900
00902 inline bool ConstIterator2D::operator!=(const ConstIterator2D& other) const
00903 {
00904 if ( msv.mdata + msv.mrange.mstart !=
00905 other.msv.mdata + other.msv.mrange.mstart )
00906 return true;
00907 else
00908 return false;
00909 }
00910
00913 inline const ConstVectorView* ConstIterator2D::operator->() const
00914 {
00915 return &msv;
00916 }
00917
00919 inline const ConstVectorView& ConstIterator2D::operator*() const
00920 {
00921 return msv;
00922 }
00923
00924
00925
00926
00927
00937 inline Index ConstVectorView::nelem() const
00938 {
00939 return mrange.mextent;
00940 }
00941
00943 inline Numeric ConstVectorView::sum() const
00944 {
00945 Numeric s=0;
00946 ConstIterator1D i = begin();
00947 const ConstIterator1D e = end();
00948
00949 for ( ; i!=e; ++i )
00950 s += *i;
00951
00952 return s;
00953 }
00954
00956 inline Numeric ConstVectorView::operator[](Index n) const
00957 {
00958
00959 assert( 0<=n );
00960 assert( n<mrange.mextent );
00961 return *( mdata +
00962 mrange.mstart +
00963 n*mrange.mstride );
00964 }
00965
00969 inline ConstVectorView ConstVectorView::operator[](const Range& r) const
00970 {
00971 return ConstVectorView(mdata, mrange, r);
00972 }
00973
00975 inline ConstIterator1D ConstVectorView::begin() const
00976 {
00977 return ConstIterator1D(mdata+mrange.mstart, mrange.mstride);
00978 }
00979
00981 inline ConstIterator1D ConstVectorView::end() const
00982 {
00983 return ConstIterator1D( mdata +
00984 mrange.mstart +
00985 (mrange.mextent)*mrange.mstride,
00986 mrange.mstride );
00987 }
00988
00990 inline ConstVectorView::operator ConstMatrixView() const
00991 {
00992 return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
00993 }
00994
00997 inline ConstVectorView::ConstVectorView() :
00998 mrange(0,0), mdata(NULL)
00999 {
01000
01001 }
01002
01005 inline ConstVectorView::ConstVectorView(Numeric *data,
01006 const Range& range) :
01007 mrange(range),
01008 mdata(data)
01009 {
01010
01011 }
01012
01023 inline ConstVectorView::ConstVectorView(Numeric *data,
01024 const Range& p,
01025 const Range& n) :
01026 mrange(p,n),
01027 mdata(data)
01028 {
01029
01030 }
01031
01035 inline std::ostream& operator<<(std::ostream& os, const ConstVectorView& v)
01036 {
01037 ConstIterator1D i=v.begin();
01038 const ConstIterator1D end=v.end();
01039
01040 if ( i!=end )
01041 {
01042 os << setw(3) << *i;
01043 ++i;
01044 }
01045 for ( ; i!=end ; ++i )
01046 {
01047 os << "\n" << setw(3) << *i;
01048 }
01049
01050 return os;
01051 }
01052
01053
01054
01055
01056
01059 inline Numeric VectorView::operator[](Index n) const
01060 {
01061 return ConstVectorView::operator[](n);
01062 }
01063
01068 inline ConstVectorView VectorView::operator[](const Range& r) const
01069 {
01070 return ConstVectorView::operator[](r);
01071 }
01072
01074 inline Numeric& VectorView::operator[](Index n)
01075 {
01076
01077 assert( 0<=n );
01078 assert( n<mrange.mextent );
01079 return *( mdata +
01080 mrange.mstart +
01081 n*mrange.mstride );
01082 }
01083
01087 inline VectorView VectorView::operator[](const Range& r)
01088 {
01089 return VectorView(mdata, mrange, r);
01090 }
01091
01095 inline ConstIterator1D VectorView::begin() const
01096 {
01097 return ConstVectorView::begin();
01098 }
01099
01103 inline ConstIterator1D VectorView::end() const
01104 {
01105 return ConstVectorView::end();
01106 }
01107
01109 inline Iterator1D VectorView::begin()
01110 {
01111 return Iterator1D(mdata+mrange.mstart, mrange.mstride);
01112 }
01113
01115 inline Iterator1D VectorView::end()
01116 {
01117 return Iterator1D( mdata +
01118 mrange.mstart +
01119 (mrange.mextent)*mrange.mstride,
01120 mrange.mstride );
01121 }
01122
01127 inline VectorView VectorView::operator=(const ConstVectorView& v)
01128 {
01129
01130
01131
01132 assert(mrange.mextent==v.mrange.mextent);
01133
01134 copy( v.begin(), v.end(), begin() );
01135
01136 return *this;
01137 }
01138
01144 inline VectorView VectorView::operator=(const VectorView& v)
01145 {
01146
01147
01148
01149 assert(mrange.mextent==v.mrange.mextent);
01150
01151 copy( v.begin(), v.end(), begin() );
01152
01153 return *this;
01154 }
01155
01158 inline VectorView VectorView::operator=(const Vector& v)
01159 {
01160
01161
01162
01163 assert(mrange.mextent==v.mrange.mextent);
01164
01165 copy( v.begin(), v.end(), begin() );
01166
01167 return *this;
01168 }
01169
01172 inline VectorView VectorView::operator=(Numeric x)
01173 {
01174 copy( x, begin(), end() );
01175 return *this;
01176 }
01177
01179 inline VectorView VectorView::operator*=(Numeric x)
01180 {
01181 const Iterator1D e=end();
01182 for ( Iterator1D i=begin(); i!=e ; ++i )
01183 *i *= x;
01184 return *this;
01185 }
01186
01188 inline VectorView VectorView::operator/=(Numeric x)
01189 {
01190 const Iterator1D e=end();
01191 for ( Iterator1D i=begin(); i!=e ; ++i )
01192 *i /= x;
01193 return *this;
01194 }
01195
01197 inline VectorView VectorView::operator+=(Numeric x)
01198 {
01199 const Iterator1D e=end();
01200 for ( Iterator1D i=begin(); i!=e ; ++i )
01201 *i += x;
01202 return *this;
01203 }
01204
01206 inline VectorView VectorView::operator-=(Numeric x)
01207 {
01208 const Iterator1D e=end();
01209 for ( Iterator1D i=begin(); i!=e ; ++i )
01210 *i -= x;
01211 return *this;
01212 }
01213
01215 inline VectorView VectorView::operator*=(const ConstVectorView& x)
01216 {
01217 assert( nelem()==x.nelem() );
01218
01219 ConstIterator1D s=x.begin();
01220
01221 Iterator1D i=begin();
01222 const Iterator1D e=end();
01223
01224 for ( ; i!=e ; ++i,++s )
01225 *i *= *s;
01226 return *this;
01227 }
01228
01230 inline VectorView VectorView::operator/=(const ConstVectorView& x)
01231 {
01232 assert( nelem()==x.nelem() );
01233
01234 ConstIterator1D s=x.begin();
01235
01236 Iterator1D i=begin();
01237 const Iterator1D e=end();
01238
01239 for ( ; i!=e ; ++i,++s )
01240 *i /= *s;
01241 return *this;
01242 }
01243
01245 inline VectorView VectorView::operator+=(const ConstVectorView& x)
01246 {
01247 assert( nelem()==x.nelem() );
01248
01249 ConstIterator1D s=x.begin();
01250
01251 Iterator1D i=begin();
01252 const Iterator1D e=end();
01253
01254 for ( ; i!=e ; ++i,++s )
01255 *i += *s;
01256 return *this;
01257 }
01258
01260 inline VectorView VectorView::operator-=(const ConstVectorView& x)
01261 {
01262 assert( nelem()==x.nelem() );
01263
01264 ConstIterator1D s=x.begin();
01265
01266 Iterator1D i=begin();
01267 const Iterator1D e=end();
01268
01269 for ( ; i!=e ; ++i,++s )
01270 *i -= *s;
01271 return *this;
01272 }
01273
01275 inline VectorView::operator MatrixView()
01276 {
01277 return MatrixView(mdata,mrange,Range(mrange.mstart,1));
01278 }
01279
01282 inline VectorView::VectorView() :
01283 ConstVectorView()
01284 {
01285
01286 }
01287
01290 inline VectorView::VectorView(Numeric *data,
01291 const Range& range) :
01292 ConstVectorView(data,range)
01293 {
01294
01295 }
01296
01307 inline VectorView::VectorView(Numeric *data,
01308 const Range& p,
01309 const Range& n) :
01310 ConstVectorView(data,p,n)
01311 {
01312
01313 }
01314
01319 inline void copy(ConstIterator1D origin,
01320 const ConstIterator1D& end,
01321 Iterator1D target)
01322 {
01323 for ( ; origin!=end ; ++origin,++target )
01324 *target = *origin;
01325 }
01326
01328 inline void copy(Numeric x,
01329 Iterator1D target,
01330 const Iterator1D& end)
01331 {
01332 for ( ; target!=end ; ++target )
01333 *target = x;
01334 }
01335
01336
01337
01338
01339
01341 inline Vector::Vector()
01342 {
01343
01344 }
01345
01347 inline Vector::Vector(Index n) :
01348 VectorView( new Numeric[n],
01349 Range(0,n))
01350 {
01351
01352 }
01353
01355 inline Vector::Vector(Index n, Numeric fill) :
01356 VectorView( new Numeric[n],
01357 Range(0,n))
01358 {
01359
01360
01361 for ( Numeric *x=mdata; x<mdata+n; ++x )
01362 *x = fill;
01363 }
01364
01373 inline Vector::Vector(Numeric start, Index extent, Numeric stride) :
01374 VectorView( new Numeric[extent],
01375 Range(0,extent))
01376 {
01377
01378 Numeric x = start;
01379 Iterator1D i=begin();
01380 const Iterator1D e=end();
01381 for ( ; i!=e; ++i )
01382 {
01383 *i = x;
01384 x += stride;
01385 }
01386 }
01387
01393 inline Vector::Vector(const ConstVectorView& v) :
01394 VectorView( new Numeric[v.nelem()],
01395 Range(0,v.nelem()))
01396 {
01397 copy(v.begin(),v.end(),begin());
01398 }
01399
01402 inline Vector::Vector(const Vector& v) :
01403 VectorView( new Numeric[v.nelem()],
01404 Range(0,v.nelem()))
01405 {
01406 copy(v.begin(),v.end(),begin());
01407 }
01408
01413 inline Vector& Vector::operator=(const Vector& v)
01414 {
01415
01416
01417
01418 assert(mrange.mextent==v.mrange.mextent);
01419 copy( v.begin(), v.end(), begin() );
01420 return *this;
01421 }
01422
01439 inline Vector& Vector::operator=(const Array<Numeric>& x)
01440 {
01441 VectorView::operator=(x);
01442 return *this;
01443 }
01444
01447 inline Vector& Vector::operator=(Numeric x)
01448 {
01449 VectorView::operator=(x);
01450 return *this;
01451 }
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01467 inline void Vector::resize(Index n)
01468 {
01469 if ( mrange.mextent != n )
01470 {
01471 delete mdata;
01472 mdata = new Numeric[n];
01473 mrange.mstart = 0;
01474 mrange.mextent = n;
01475 mrange.mstride = 1;
01476 }
01477 }
01478
01481 inline Vector::~Vector()
01482 {
01483 delete [] mdata;
01484 }
01485
01486
01487
01488
01489
01491 inline Index ConstMatrixView::nrows() const
01492 {
01493 return mrr.mextent;
01494 }
01495
01497 inline Index ConstMatrixView::ncols() const
01498 {
01499 return mcr.mextent;
01500 }
01501
01503 inline Numeric ConstMatrixView::operator()(Index r, Index c) const
01504 {
01505
01506 assert( 0<=r );
01507 assert( 0<=c );
01508 assert( r<mrr.mextent );
01509 assert( c<mcr.mextent );
01510
01511 return *( mdata +
01512 mrr.mstart +
01513 r*mrr.mstride +
01514 mcr.mstart +
01515 c*mcr.mstride );
01516 }
01517
01521 inline ConstMatrixView ConstMatrixView::operator()(const Range& r,
01522 const Range& c) const
01523 {
01524 return ConstMatrixView(mdata, mrr, mcr, r, c);
01525 }
01526
01532 inline ConstVectorView ConstMatrixView::operator()(const Range& r, Index c) const
01533 {
01534
01535 assert( 0 <= c );
01536 assert( c < mcr.mextent );
01537
01538 return ConstVectorView(mdata + mcr.mstart + c*mcr.mstride,
01539 mrr, r);
01540 }
01541
01547 inline ConstVectorView ConstMatrixView::operator()(Index r, const Range& c) const
01548 {
01549
01550 assert( 0 <= r );
01551 assert( r < mrr.mextent );
01552
01553 return ConstVectorView(mdata + mrr.mstart + r*mrr.mstride,
01554 mcr, c);
01555 }
01556
01558 inline ConstIterator2D ConstMatrixView::begin() const
01559 {
01560 return ConstIterator2D(ConstVectorView(mdata+mrr.mstart,
01561 mcr),
01562 mrr.mstride);
01563 }
01564
01566 inline ConstIterator2D ConstMatrixView::end() const
01567 {
01568 return ConstIterator2D( ConstVectorView(mdata + mrr.mstart + (mrr.mextent)*mrr.mstride,
01569 mcr),
01570 mrr.mstride );
01571 }
01572
01575 inline ConstMatrixView::ConstMatrixView() :
01576 mrr(0,0,1), mcr(0,0,1), mdata(NULL)
01577 {
01578
01579 }
01580
01584 inline ConstMatrixView::ConstMatrixView(Numeric *data,
01585 const Range& rr,
01586 const Range& cr) :
01587 mrr(rr),
01588 mcr(cr),
01589 mdata(data)
01590 {
01591
01592 }
01593
01605 inline ConstMatrixView::ConstMatrixView(Numeric *data,
01606 const Range& pr, const Range& pc,
01607 const Range& nr, const Range& nc) :
01608 mrr(pr,nr),
01609 mcr(pc,nc),
01610 mdata(data)
01611 {
01612
01613 }
01614
01622 inline std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v)
01623 {
01624
01625 ConstIterator2D ir=v.begin();
01626 const ConstIterator2D end_row=v.end();
01627
01628 if ( ir!=end_row )
01629 {
01630 ConstIterator1D ic = ir->begin();
01631 const ConstIterator1D end_col = ir->end();
01632
01633 if ( ic!=end_col )
01634 {
01635 os << setw(3) << *ic;
01636 ++ic;
01637 }
01638 for ( ; ic!=end_col ; ++ic )
01639 {
01640 os << " " << setw(3) << *ic;
01641 }
01642 ++ir;
01643 }
01644 for ( ; ir!=end_row ; ++ir )
01645 {
01646 ConstIterator1D ic = ir->begin();
01647 const ConstIterator1D end_col = ir->end();
01648
01649 os << "\n";
01650 if ( ic!=end_col )
01651 {
01652 os << setw(3) << *ic;
01653 ++ic;
01654 }
01655 for ( ; ic!=end_col ; ++ic )
01656 {
01657 os << " " << setw(3) << *ic;
01658 }
01659 }
01660
01661 return os;
01662 }
01663
01664
01665
01666
01667
01670 inline Numeric MatrixView::operator()(Index r, Index c) const
01671 {
01672 return ConstMatrixView::operator()(r,c);
01673 }
01674
01679 inline ConstMatrixView MatrixView::operator()(const Range& r, const Range& c) const
01680 {
01681 return ConstMatrixView::operator()(r,c);
01682 }
01683
01690 inline ConstVectorView MatrixView::operator()(const Range& r, Index c) const
01691 {
01692 return ConstMatrixView::operator()(r,c);
01693 }
01694
01701 inline ConstVectorView MatrixView::operator()(Index r, const Range& c) const
01702 {
01703 return ConstMatrixView::operator()(r,c);
01704 }
01705
01707 inline Numeric& MatrixView::operator()(Index r, Index c)
01708 {
01709
01710 assert( 0<=r );
01711 assert( 0<=c );
01712 assert( r<mrr.mextent );
01713 assert( c<mcr.mextent );
01714
01715 return *( mdata +
01716 mrr.mstart +
01717 r*mrr.mstride +
01718 mcr.mstart +
01719 c*mcr.mstride );
01720 }
01721
01725 inline MatrixView MatrixView::operator()(const Range& r, const Range& c)
01726 {
01727 return MatrixView(mdata, mrr, mcr, r, c);
01728 }
01729
01734 inline VectorView MatrixView::operator()(const Range& r, Index c)
01735 {
01736
01737 assert( 0 <= c );
01738 assert( c < mcr.mextent );
01739
01740 return VectorView(mdata + mcr.mstart + c*mcr.mstride,
01741 mrr, r);
01742 }
01743
01748 inline VectorView MatrixView::operator()(Index r, const Range& c)
01749 {
01750
01751 assert( 0 <= r );
01752 assert( r < mrr.mextent );
01753
01754 return VectorView(mdata + mrr.mstart + r*mrr.mstride,
01755 mcr, c);
01756 }
01757
01760 inline ConstIterator2D MatrixView::begin() const
01761 {
01762 return ConstMatrixView::begin();
01763 }
01764
01766 inline ConstIterator2D MatrixView::end() const
01767 {
01768 return ConstMatrixView::end();
01769 }
01770
01772 inline Iterator2D MatrixView::begin()
01773 {
01774 return Iterator2D(VectorView(mdata+mrr.mstart, mcr),
01775 mrr.mstride);
01776 }
01777
01779 inline Iterator2D MatrixView::end()
01780 {
01781 return Iterator2D( VectorView(mdata + mrr.mstart + (mrr.mextent)*mrr.mstride,
01782 mcr),
01783 mrr.mstride );
01784 }
01785
01790 inline MatrixView& MatrixView::operator=(const ConstMatrixView& m)
01791 {
01792
01793 assert(mrr.mextent==m.mrr.mextent);
01794 assert(mcr.mextent==m.mcr.mextent);
01795
01796 copy( m.begin(), m.end(), begin() );
01797 return *this;
01798 }
01799
01805 inline MatrixView& MatrixView::operator=(const MatrixView& m)
01806 {
01807
01808 assert(mrr.mextent==m.mrr.mextent);
01809 assert(mcr.mextent==m.mcr.mextent);
01810
01811 copy( m.begin(), m.end(), begin() );
01812 return *this;
01813 }
01814
01818 inline MatrixView& MatrixView::operator=(const Matrix& m)
01819 {
01820
01821 assert(mrr.mextent==m.mrr.mextent);
01822 assert(mcr.mextent==m.mcr.mextent);
01823
01824 copy( m.begin(), m.end(), begin() );
01825 return *this;
01826 }
01827
01832 inline MatrixView& MatrixView::operator=(const ConstVectorView& v)
01833 {
01834
01835 assert( mrr.mextent==v.nelem() );
01836 assert( mcr.mextent==1 );
01837
01838 ConstMatrixView dummy(v);
01839 copy( dummy.begin(), dummy.end(), begin() );
01840 return *this;
01841 }
01842
01845 inline MatrixView& MatrixView::operator=(Numeric x)
01846 {
01847 copy( x, begin(), end() );
01848 return *this;
01849 }
01850
01852 inline MatrixView& MatrixView::operator*=(Numeric x)
01853 {
01854 const Iterator2D er=end();
01855 for ( Iterator2D r=begin(); r!=er ; ++r )
01856 {
01857 const Iterator1D ec = r->end();
01858 for ( Iterator1D c = r->begin(); c!=ec ; ++c )
01859 *c *= x;
01860 }
01861 return *this;
01862 }
01863
01865 inline MatrixView& MatrixView::operator/=(Numeric x)
01866 {
01867 const Iterator2D er=end();
01868 for ( Iterator2D r=begin(); r!=er ; ++r )
01869 {
01870 const Iterator1D ec = r->end();
01871 for ( Iterator1D c = r->begin(); c!=ec ; ++c )
01872 *c /= x;
01873 }
01874 return *this;
01875 }
01876
01878 inline MatrixView& MatrixView::operator+=(Numeric x)
01879 {
01880 const Iterator2D er=end();
01881 for ( Iterator2D r=begin(); r!=er ; ++r )
01882 {
01883 const Iterator1D ec = r->end();
01884 for ( Iterator1D c = r->begin(); c!=ec ; ++c )
01885 *c += x;
01886 }
01887 return *this;
01888 }
01889
01891 inline MatrixView& MatrixView::operator-=(Numeric x)
01892 {
01893 const Iterator2D er=end();
01894 for ( Iterator2D r=begin(); r!=er ; ++r )
01895 {
01896 const Iterator1D ec = r->end();
01897 for ( Iterator1D c = r->begin(); c!=ec ; ++c )
01898 *c -= x;
01899 }
01900 return *this;
01901 }
01902
01904 inline MatrixView& MatrixView::operator*=(const ConstMatrixView& x)
01905 {
01906 assert(nrows()==x.nrows());
01907 assert(ncols()==x.ncols());
01908 ConstIterator2D sr = x.begin();
01909 Iterator2D r = begin();
01910 const Iterator2D er = end();
01911 for ( ; r!=er ; ++r,++sr )
01912 {
01913 ConstIterator1D sc = sr->begin();
01914 Iterator1D c = r->begin();
01915 const Iterator1D ec = r->end();
01916 for ( ; c!=ec ; ++c,++sc )
01917 *c *= *sc;
01918 }
01919 return *this;
01920 }
01921
01923 inline MatrixView& MatrixView::operator/=(const ConstMatrixView& x)
01924 {
01925 assert(nrows()==x.nrows());
01926 assert(ncols()==x.ncols());
01927 ConstIterator2D sr = x.begin();
01928 Iterator2D r = begin();
01929 const Iterator2D er = end();
01930 for ( ; r!=er ; ++r,++sr )
01931 {
01932 ConstIterator1D sc = sr->begin();
01933 Iterator1D c = r->begin();
01934 const Iterator1D ec = r->end();
01935 for ( ; c!=ec ; ++c,++sc )
01936 *c /= *sc;
01937 }
01938 return *this;
01939 }
01940
01942 inline MatrixView& MatrixView::operator+=(const ConstMatrixView& x)
01943 {
01944 assert(nrows()==x.nrows());
01945 assert(ncols()==x.ncols());
01946 ConstIterator2D sr = x.begin();
01947 Iterator2D r = begin();
01948 const Iterator2D er = end();
01949 for ( ; r!=er ; ++r,++sr )
01950 {
01951 ConstIterator1D sc = sr->begin();
01952 Iterator1D c = r->begin();
01953 const Iterator1D ec = r->end();
01954 for ( ; c!=ec ; ++c,++sc )
01955 *c += *sc;
01956 }
01957 return *this;
01958 }
01959
01961 inline MatrixView& MatrixView::operator-=(const ConstMatrixView& x)
01962 {
01963 assert(nrows()==x.nrows());
01964 assert(ncols()==x.ncols());
01965 ConstIterator2D sr = x.begin();
01966 Iterator2D r = begin();
01967 const Iterator2D er = end();
01968 for ( ; r!=er ; ++r,++sr )
01969 {
01970 ConstIterator1D sc = sr->begin();
01971 Iterator1D c = r->begin();
01972 const Iterator1D ec = r->end();
01973 for ( ; c!=ec ; ++c,++sc )
01974 *c -= *sc;
01975 }
01976 return *this;
01977 }
01978
01980 inline MatrixView& MatrixView::operator*=(const ConstVectorView& x)
01981 {
01982 assert(nrows()==x.nelem());
01983 assert(ncols()==1);
01984 ConstIterator1D sc = x.begin();
01985 Iterator2D r = begin();
01986 const Iterator2D er = end();
01987 for ( ; r!=er ; ++r,++sc )
01988 {
01989 Iterator1D c = r->begin();
01990 *c *= *sc;
01991 }
01992 return *this;
01993 }
01994
01996 inline MatrixView& MatrixView::operator/=(const ConstVectorView& x)
01997 {
01998 assert(nrows()==x.nelem());
01999 assert(ncols()==1);
02000 ConstIterator1D sc = x.begin();
02001 Iterator2D r = begin();
02002 const Iterator2D er = end();
02003 for ( ; r!=er ; ++r,++sc )
02004 {
02005 Iterator1D c = r->begin();
02006 *c /= *sc;
02007 }
02008 return *this;
02009 }
02010
02012 inline MatrixView& MatrixView::operator+=(const ConstVectorView& x)
02013 {
02014 assert(nrows()==x.nelem());
02015 assert(ncols()==1);
02016 ConstIterator1D sc = x.begin();
02017 Iterator2D r = begin();
02018 const Iterator2D er = end();
02019 for ( ; r!=er ; ++r,++sc )
02020 {
02021 Iterator1D c = r->begin();
02022 *c += *sc;
02023 }
02024 return *this;
02025 }
02026
02028 inline MatrixView& MatrixView::operator-=(const ConstVectorView& x)
02029 {
02030 assert(nrows()==x.nelem());
02031 assert(ncols()==1);
02032 ConstIterator1D sc = x.begin();
02033 Iterator2D r = begin();
02034 const Iterator2D er = end();
02035 for ( ; r!=er ; ++r,++sc )
02036 {
02037 Iterator1D c = r->begin();
02038 *c -= *sc;
02039 }
02040 return *this;
02041 }
02042
02045 inline MatrixView::MatrixView() :
02046 ConstMatrixView()
02047 {
02048
02049 }
02050
02054 inline MatrixView::MatrixView(Numeric *data,
02055 const Range& rr,
02056 const Range& cr) :
02057 ConstMatrixView(data, rr, cr)
02058 {
02059
02060 }
02061
02073 inline MatrixView::MatrixView(Numeric *data,
02074 const Range& pr, const Range& pc,
02075 const Range& nr, const Range& nc) :
02076 ConstMatrixView(data,pr,pc,nr,nc)
02077 {
02078
02079 }
02080
02090 inline void copy(ConstIterator2D origin,
02091 const ConstIterator2D& end,
02092 Iterator2D target)
02093 {
02094 for ( ; origin!=end ; ++origin,++target )
02095 {
02096 ConstIterator1D o = origin->begin();
02097 const ConstIterator1D e = origin->end();
02098 Iterator1D t = target->begin();
02099 for ( ; o!=e ; ++o,++t )
02100 *t = *o;
02101 }
02102 }
02103
02105 inline void copy(Numeric x,
02106 Iterator2D target,
02107 const Iterator2D& end)
02108 {
02109 for ( ; target!=end ; ++target )
02110 {
02111 Iterator1D t = target->begin();
02112 const Iterator1D e = target->end();
02113 for ( ; t!=e ; ++t )
02114 *t = x;
02115 }
02116 }
02117
02118
02119
02120
02121
02123 inline Matrix::Matrix() :
02124 MatrixView::MatrixView()
02125 {
02126
02127
02128
02129
02130 }
02131
02134 inline Matrix::Matrix(Index r, Index c) :
02135 MatrixView( new Numeric[r*c],
02136 Range(0,r,c),
02137 Range(0,c))
02138 {
02139
02140 }
02141
02143 inline Matrix::Matrix(Index r, Index c, Numeric fill) :
02144 MatrixView( new Numeric[r*c],
02145 Range(0,r,c),
02146 Range(0,c))
02147 {
02148
02149
02150 for ( Numeric *x=mdata; x<mdata+r*c; ++x )
02151 *x = fill;
02152 }
02153
02156 inline Matrix::Matrix(const ConstMatrixView& m) :
02157 MatrixView( new Numeric[m.nrows()*m.ncols()],
02158 Range( 0, m.nrows(), m.ncols() ),
02159 Range( 0, m.ncols() ) )
02160 {
02161 copy(m.begin(),m.end(),begin());
02162 }
02163
02166 inline Matrix::Matrix(const Matrix& m) :
02167 MatrixView( new Numeric[m.nrows()*m.ncols()],
02168 Range( 0, m.nrows(), m.ncols() ),
02169 Range( 0, m.ncols() ) )
02170 {
02171
02172
02173
02174
02175 copy(m.begin(),m.end(),begin());
02176 }
02177
02188 inline Matrix& Matrix::operator=(const Matrix& m)
02189 {
02190
02191
02192
02193
02194
02195 if ( 0 == mrr.mextent )
02196 {
02197
02198 resize( m.mrr.mextent, m.mcr.mextent );
02199 }
02200 else
02201 {
02202
02203 assert( mrr.mextent==m.mrr.mextent );
02204 assert( mcr.mextent==m.mcr.mextent );
02205 }
02206
02207 copy( m.begin(), m.end(), begin() );
02208 return *this;
02209 }
02210
02213 inline Matrix& Matrix::operator=(Numeric x)
02214 {
02215 copy( x, begin(), end() );
02216 return *this;
02217 }
02218
02223 inline Matrix& Matrix::operator=(const ConstVectorView& v)
02224 {
02225
02226 assert( mrr.mextent==v.nelem() );
02227 assert( mcr.mextent==1 );
02228
02229 ConstMatrixView dummy(v);
02230 copy( dummy.begin(), dummy.end(), begin() );
02231 return *this;
02232 }
02233
02237 inline void Matrix::resize(Index r, Index c)
02238 {
02239 assert( 0<=r );
02240 assert( 0<=c );
02241
02242 if ( mrr.mextent!=r || mcr.mextent!=c )
02243 {
02244 delete mdata;
02245 mdata = new Numeric[r*c];
02246
02247 mrr.mstart = 0;
02248 mrr.mextent = r;
02249 mrr.mstride = c;
02250
02251 mcr.mstart = 0;
02252 mcr.mextent = c;
02253 mcr.mstride = 1;
02254 }
02255 }
02256
02259 inline Matrix::~Matrix()
02260 {
02261
02262
02263 delete [] mdata;
02264 }
02265
02266
02267
02268
02270 inline Numeric operator*(const ConstVectorView& a, const ConstVectorView& b)
02271 {
02272
02273 assert( a.nelem() == b.nelem() );
02274
02275 const ConstIterator1D ae = a.end();
02276 ConstIterator1D ai = a.begin();
02277 ConstIterator1D bi = b.begin();
02278
02279 Numeric res = 0;
02280 for ( ; ai!=ae ; ++ai, ++bi )
02281 res += (*ai) * (*bi);
02282
02283 return res;
02284 }
02285
02291 inline void mult( VectorView y,
02292 const ConstMatrixView& M,
02293 const ConstVectorView& x )
02294 {
02295
02296 assert( y.nelem() == M.nrows() );
02297 assert( M.ncols() == x.nelem() );
02298
02299
02300 const Iterator1D ye = y.end();
02301 Iterator1D yi = y.begin();
02302
02303
02304 ConstIterator2D mi = M.begin();
02305
02306
02307 const ConstIterator1D xs = x.begin();
02308 const ConstIterator1D xe = x.end();
02309
02310
02311 for ( ; yi!=ye ; ++yi, ++mi )
02312 {
02313
02314 Numeric dummy = 0;
02315
02316
02317 ConstIterator1D ri = mi->begin();
02318
02319
02320 ConstIterator1D xi = xs;
02321
02322 for ( ; xi!=xe ; ++xi, ++ri )
02323 {
02324 dummy += (*ri) * (*xi);
02325 }
02326
02327 *yi = dummy;
02328 }
02329 }
02330
02336 inline void mult( MatrixView A,
02337 const ConstMatrixView& B,
02338 const ConstMatrixView& C )
02339 {
02340
02341 assert( A.nrows() == B.nrows() );
02342 assert( A.ncols() == C.ncols() );
02343 assert( B.ncols() == C.nrows() );
02344
02345
02346
02347 ConstMatrixView CT = transpose(C);
02348
02349 const Iterator2D ae = A.end();
02350 Iterator2D ai = A.begin();
02351 ConstIterator2D bi = B.begin();
02352
02353
02354 for ( ; ai!=ae ; ++ai, ++bi )
02355 {
02356 const Iterator1D ace = ai->end();
02357 Iterator1D aci = ai->begin();
02358 ConstIterator2D cti = CT.begin();
02359
02360
02361
02362
02363 for ( ; aci!=ace ; ++aci, ++cti )
02364 {
02365
02366
02367 *aci = (*bi) * (*cti);
02368 }
02369 }
02370 }
02371
02373 inline ConstMatrixView transpose(ConstMatrixView m)
02374 {
02375 return ConstMatrixView(m.mdata, m.mcr, m.mrr);
02376 }
02377
02380 inline MatrixView transpose(MatrixView m)
02381 {
02382 return MatrixView(m.mdata, m.mcr, m.mrr);
02383 }
02384
02405 inline void transform( VectorView y,
02406 double (&my_func)(double),
02407 ConstVectorView x )
02408 {
02409
02410 assert( y.nelem()==x.nelem() );
02411
02412 const ConstIterator1D xe = x.end();
02413 ConstIterator1D xi = x.begin();
02414 Iterator1D yi = y.begin();
02415 for ( ; xi!=xe; ++xi, ++yi )
02416 *yi = my_func(*xi);
02417 }
02418
02437 inline void transform( MatrixView y,
02438 double (&my_func)(double),
02439 ConstMatrixView x )
02440 {
02441
02442 assert( y.nrows()==x.nrows() );
02443 assert( y.ncols()==x.ncols() );
02444
02445 const ConstIterator2D rxe = x.end();
02446 ConstIterator2D rx = x.begin();
02447 Iterator2D ry = y.begin();
02448 for ( ; rx!=rxe; ++rx, ++ry )
02449 {
02450 const ConstIterator1D cxe = rx->end();
02451 ConstIterator1D cx = rx->begin();
02452 Iterator1D cy = ry->begin();
02453 for ( ; cx!=cxe; ++cx, ++cy )
02454 *cy = my_func(*cx);
02455 }
02456 }
02457
02459 inline Numeric max(const ConstVectorView& x)
02460 {
02461
02462 Numeric max = x[0];
02463
02464 const ConstIterator1D xe = x.end();
02465 ConstIterator1D xi = x.begin();
02466
02467 for ( ; xi!=xe ; ++xi )
02468 {
02469 if ( *xi > max )
02470 max = *xi;
02471 }
02472
02473 return max;
02474 }
02475
02477 inline Numeric max(const ConstMatrixView& x)
02478 {
02479
02480 Numeric max = x(0,0);
02481
02482 const ConstIterator2D rxe = x.end();
02483 ConstIterator2D rx = x.begin();
02484
02485 for ( ; rx!=rxe ; ++rx )
02486 {
02487 const ConstIterator1D cxe = rx->end();
02488 ConstIterator1D cx = rx->begin();
02489
02490 for ( ; cx!=cxe ; ++cx )
02491 if ( *cx > max )
02492 max = *cx;
02493 }
02494
02495 return max;
02496 }
02497
02499 inline Numeric min(const ConstVectorView& x)
02500 {
02501
02502 Numeric min = x[0];
02503
02504 const ConstIterator1D xe = x.end();
02505 ConstIterator1D xi = x.begin();
02506
02507 for ( ; xi!=xe ; ++xi )
02508 {
02509 if ( *xi < min )
02510 min = *xi;
02511 }
02512
02513 return min;
02514 }
02515
02517 inline Numeric min(const ConstMatrixView& x)
02518 {
02519
02520 Numeric min = x(0,0);
02521
02522 const ConstIterator2D rxe = x.end();
02523 ConstIterator2D rx = x.begin();
02524
02525 for ( ; rx!=rxe ; ++rx )
02526 {
02527 const ConstIterator1D cxe = rx->end();
02528 ConstIterator1D cx = rx->begin();
02529
02530 for ( ; cx!=cxe ; ++cx )
02531 if ( *cx < min )
02532 min = *cx;
02533 }
02534
02535 return min;
02536 }
02537
02538
02539 #endif // matpackI_h