more/num/matrix.h Source File

00001 //  Copyright (C) 1998--2002  Petter Urkedal
00002 //
00003 //  This file is free software; you can redistribute it and/or modify
00004 //  it under the terms of the GNU General Public License as published by
00005 //  the Free Software Foundation; either version 2 of the License, or
00006 //  (at your option) any later version.
00007 //
00008 //  This file is distributed in the hope that it will be useful,
00009 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011 //  GNU General Public License for more details.
00012 //
00013 //  You should have received a copy of the GNU General Public License
00014 //  along with this program; if not, write to the Free Software
00015 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00016 //
00017 //  As a special exception, you may use this file as part of a free
00018 //  software library without restriction.  Specifically, if other files
00019 //  instantiate templates or use macros or inline functions from this
00020 //  file, or you compile this file and link it with other files to
00021 //  produce an executable, this file does not by itself cause the
00022 //  resulting executable to be covered by the GNU General Public
00023 //  License.  This exception does not however invalidate any other
00024 //  reasons why the executable file might be covered by the GNU General
00025 //  Public License.
00026 
00027 //  $Id: matrix.h,v 1.2 2002/09/01 12:02:10 petter_urkedal Exp $
00028 
00029 
00030 #ifndef MORE_NUM_MATRIX_H
00031 #define MORE_NUM_MATRIX_H
00032 
00033 ///\if unmaintained
00034 #include <complex>
00035 #include <iterator>
00036 #include <more/bits/matrixtypes.h>
00037 #include <more/bits/mathbits.h>
00038 #include <more/gen/utility.h>
00039 #include <more/meta.h>
00040 
00041 
00042 namespace more {
00043 namespace num {
00044 
00045   //  ---   M a t r i x   V i e w s   ---
00046   //
00047   //  These are utilities when defining matrices.  They should not be
00048   //  used directly, but accessed through the typedefs of matrix
00049   //  classes.
00050 
00051   template <class T, class Signature = matrix_s>
00052     class rectangular_view
00053     {
00054         int size0;
00055         T* array;
00056     public:
00057         rectangular_view(T& x, int n) : array(&x), size0(n) {}
00058         T& operator()(int i, int j) { return array[i+j*size0]; }
00059         T operator()(int i, int j) const { return array[i+j*size0]; }
00060     };
00061 
00062   template <class T>
00063     class vector_view
00064     {
00065         int stride;
00066         T* array;
00067     public:
00068         vector_view(T& x, int n) : array(&x), stride(n) {}
00069         T& operator()(int i) { return array[i*stride]; }
00070         T operator()(int i) const { return array[i*stride]; }
00071     };
00072 
00073 
00074 
00075     //  ---   A   C o n v e n i e n t   I t e r a t o r   C l a s s   ---
00076 
00077   template <class T, class Difference = ptrdiff_t>
00078     struct rank1_iterator
00079         : std::iterator<std::random_access_iterator_tag, T, Difference>
00080     {
00081         typedef std::iterator<std::random_access_iterator_tag, T, Difference>
00082                 base;
00083         typedef rank1_iterator self;
00084         typedef typename base::reference reference;
00085         typedef typename base::pointer pointer;
00086         typedef typename base::difference_type difference_type;
00087         rank1_iterator() {}
00088         rank1_iterator(const rank1_iterator& it)
00089             : stride(it.stride), current(it.current) {}
00090         rank1_iterator(T* ptr, int stride_)
00091             : stride(stride_), current(ptr) {}
00092         bool operator==(const self& x) const { return current==x.current; }
00093         bool operator!=(const self& x) const { return current!=x.current; }
00094         bool operator<(const self& x) const { return current<x.current; }
00095         reference operator*() { return *current; }
00096         self& operator++() { current += stride; return *this; }
00097         self operator++(int)
00098             { self tmp = *this; current += stride; return tmp; }
00099         self& operator--() { current -= stride; return *this; }
00100         self operator--(int)
00101             { self tmp = *this; current -= stride; return tmp; }
00102         self operator+(difference_type n) const
00103             { return self(current+n*stride, stride); }
00104         self operator-(difference_type n) const
00105             { return self(current-n*stride, stride); }
00106         difference_type operator-(const self& rhs) const
00107             { return (current - rhs.current)/stride; }
00108         self& operator+=(difference_type n)
00109             { current += n*stride; return *this; }
00110         self& operator-=(difference_type n)
00111             { current -= n*stride; return *this; }
00112         reference operator[](difference_type n)
00113             {return current[n*stride];}
00114       protected:
00115         int stride;
00116         T *current;
00117     };
00118 
00119   template<typename T>
00120     struct rank2_proxy {
00121 //      template<typename U> friend class rank2_proxy<U>;
00122         typedef rank1_iterator<T> iterator;
00123         rank2_proxy() {}
00124         template<typename U>
00125           rank2_proxy(rank2_proxy<U> const& px)
00126               : current(px.current), stride(px.stride),
00127                 substride(px.substride), subsize(px.subsize) {}
00128         rank2_proxy(T *cur, int stride_, int substride_, int subsize_)
00129             : current(cur), stride(stride_), substride(substride_),
00130               subsize(subsize_) {}
00131         rank1_iterator<T> begin()
00132             { return rank1_iterator<T>(current, substride); }
00133         rank1_iterator<T> end()
00134             { return rank1_iterator<T>(current+subsize, substride); }
00135       protected:
00136         T *current;
00137         int stride, substride, subsize;
00138     };
00139 
00140   template<typename T, typename Difference = ptrdiff_t>
00141     struct rank2_iterator
00142         : private rank2_proxy<T>,
00143           std::iterator< std::random_access_iterator_tag, rank2_proxy<T>,
00144                          Difference >
00145     {
00146 //      template<typename U, typename Diff2>
00147 //        friend class rank2_iterator<U, Diff2>;
00148         typedef std::iterator< std::random_access_iterator_tag,
00149                                rank2_proxy<T>, Difference > base;
00150         typedef rank2_iterator self;
00151         typedef rank1_iterator<T> iterator;
00152         typedef typename base::reference reference;
00153         typedef typename base::pointer pointer;
00154         typedef typename base::value_type value_type;
00155         typedef typename base::difference_type difference_type;
00156         rank2_iterator() {}
00157         rank2_iterator(rank2_iterator const& it) : rank2_proxy<T>(it) {}
00158         template<typename U, typename Diff2>
00159           rank2_iterator(rank2_iterator<U, Diff2> const& it)
00160               : rank2_proxy<T>(it._priv_base_access()) {}
00161 //      rank2_iterator(rank<typename mutable_<T>::eval, Difference> const& it)
00162 //          : rank2_proxy<T>(it.ptr, it.stride, it.substride, it.subsize) {}
00163         rank2_iterator(T* ptr, int stride_, int substride_, int subsize_)
00164             : rank2_proxy<T>(ptr, stride_, substride_, subsize_) {}
00165         bool operator==(const self& x) const { return current==x.current; }
00166         bool operator!=(const self& x) const { return current!=x.current; }
00167         bool operator<(const self& x) const { return current<x.current; }
00168         self& operator++() { current += stride; return *this; }
00169         self operator++(int)
00170             { self tmp = *this; current += stride; return tmp; }
00171         self& operator--() { current -= stride; return *this; }
00172         self operator--(int)
00173             { self tmp = *this; current -= stride; return tmp; }
00174         self operator+(difference_type n) const
00175             { return self(current+n*stride, stride, substride, subsize); }
00176         self operator-(difference_type n) const
00177             { return self(current-n*stride, stride, substride, subsize); }
00178         difference_type operator-(const self& rhs) const
00179             { return (current - rhs.current)/stride; }
00180         self& operator+=(difference_type n)
00181             { current += n*stride; return *this; }
00182         self& operator-=(difference_type n)
00183             { current -= n*stride; return *this; }
00184         reference operator*() { return *this; }
00185         pointer operator->() { return this; }
00186         value_type operator[](difference_type n)
00187             { return value_type(current+n*stride, stride, substride); }
00188         rank2_proxy<T> const& _priv_base_access() const { return *this; }
00189     };
00190 
00191 
00192   template <class Matrix, class Difference = ptrdiff_t>
00193     struct indexing_rank1_iterator
00194         : std::iterator< std::random_access_iterator_tag,
00195                          typename Matrix::value_type, Difference >
00196     {
00197         typedef typename Matrix::value_type value_type;
00198         typedef std::iterator< std::random_access_iterator_tag, value_type,
00199                                Difference > base;
00200         typedef indexing_rank1_iterator self;
00201         typedef typename base::reference reference;
00202         typedef typename base::pointer pointer;
00203         typedef typename base::difference_type difference_type;
00204         typedef typename Matrix::value_proxy value_proxy;
00205         indexing_rank1_iterator() {}
00206         indexing_rank1_iterator(const indexing_rank1_iterator& it)
00207             : i(it.i), j(it.j), M(it.M) {}
00208         indexing_rank1_iterator(Matrix& M_, int i_, int j_)
00209             : i(i_), j(j_), M(&M_) {}
00210         self& operator=(const self& rhs) { i = rhs.i; j = rhs.j; M = rhs.M; }
00211         bool operator==(const self& rhs) const { return i == rhs.i; }
00212         bool operator!=(const self& rhs) const { return i != rhs.i; }
00213         bool operator<(const self& rhs) const { return i < rhs.i; }
00214         typename if_< typename is_const_<Matrix>::eval,
00215                       value_type, value_proxy >::eval
00216           operator*() { return (*M)(i, j); }
00217         self& operator++() { ++i; return *this; }
00218         self operator++(int) { self tmp = *this; ++i; return tmp; }
00219         self& operator--() { --i; return *this; }
00220         self operator--(int) { self tmp = *this; --i; return tmp; }
00221         self operator+(difference_type n) const { return self(*M, i+n, j); }
00222         friend self operator+(difference_type n, self const& it)
00223             { return self(*it.M, it.i+n, it.j); }
00224         self operator-(difference_type n) const { return self(*M, i-n, j); }
00225         difference_type operator-(const self& rhs) const
00226             { return (i - rhs.i); }
00227         self& operator+=(difference_type n) { i += n; return *this; }
00228         self& operator-=(difference_type n) { i -= n; return *this; }
00229         reference operator[](difference_type n) const { return (*M)(i+n, j); }
00230       protected:
00231         int i, j;
00232         mutable Matrix* M;
00233     };
00234 
00235   template<typename Matrix, typename Difference = ptrdiff_t>
00236     struct indexing_rank2_iterator
00237         : std::iterator< std::random_access_iterator_tag,
00238                          indexing_rank2_iterator<Matrix, Difference>,
00239                          Difference >
00240     {
00241         typedef indexing_rank1_iterator<Matrix, Difference> iterator;
00242         typedef std::iterator< std::random_access_iterator_tag,
00243                                indexing_rank2_iterator<Matrix, Difference>,
00244                                Difference > base;
00245         typedef indexing_rank2_iterator self;
00246         typedef typename base::reference reference;
00247         typedef typename base::pointer pointer;
00248         typedef typename base::difference_type difference_type;
00249         typedef typename base::value_type value_type;
00250         indexing_rank2_iterator() {}
00251         indexing_rank2_iterator(const indexing_rank2_iterator& it)
00252             : n(it.n), j(it.j), M(it.M) {}
00253         indexing_rank2_iterator(Matrix& M_, int n_, int j_)
00254             : n(n_), j(j_), M(M_) {}
00255         bool operator==(const self& rhs) const { return j == rhs.j; }
00256         bool operator!=(const self& rhs) const { return j != rhs.j; }
00257         bool operator<(const self& rhs) const { return j < rhs.j; }
00258         reference operator*() { return *this; }
00259         pointer operator->() { return this; }
00260         self& operator++() { ++j; return *this; }
00261         self operator++(int) { self tmp = *this; ++j; return tmp; }
00262         self& operator--() { --j; return *this; }
00263         self operator--(int) { self tmp = *this; --j; return tmp; }
00264         self operator+(difference_type k) const { return self(M, n, j+k); }
00265         self operator-(difference_type k) const { return self(M, n, j-k); }
00266         difference_type operator-(const self& rhs) const
00267             { return (j - rhs.j); }
00268         self& operator+=(difference_type k) { j += k; return *this; }
00269         self& operator-=(difference_type k) { j -= k; return *this; }
00270         value_type operator[](difference_type k) const
00271             { return self(M, n, j+k); }
00272         iterator begin() { return iterator(M, 0, j); }
00273         iterator end() { return iterator(M, n, j); }
00274       protected:
00275         int n, j;
00276         Matrix& M;
00277     };
00278 
00279 
00280   template <class T, class Signature = matrix_s> class matrix;
00281 
00282   template <class T, class Signature = matrix_s>
00283     class matrix_view {
00284     public:
00285         typedef typename matrix<T, Signature>::iterator iterator;
00286         typedef iterator column_iterator;
00287         typedef rank1_iterator<T> row_iterator;
00288         typedef rank1_iterator<T> diagonal_iterator;
00289     private:
00290         int n, m, r;
00291         T* array;
00292     public:
00293         matrix_view(T* ptr, int n_, int m_, int n_big) : array(ptr)
00294             { n = n_, m = m_; r = n_big; }
00295         T operator()(int i, int j) const { return array[i+r*j]; }
00296         T& operator()(int i, int j) { return array[i+r*j]; }
00297         int size(int i) { return i==0? n : i==1? m : r; }
00298         void become_mutable() {}
00299         // etc.
00300     };
00301 
00302 
00303 
00304 //  ---   S t a n d a r d   D e n s e   M a t r i c e s   ---
00305 //
00306 //  The definitions below defines a standard dense matrix, and derives
00307 //  matrices of definite properties from this.  The standard dense
00308 //  matrix is made to be Fortran 77 compatible, but this fact should
00309 //  be utilized only when extending the present library.  A user-code
00310 //  which requires fortran compatible matrices, should use
00311 //  typedef-aliases below.  The representation of the standard
00312 //  matrices may be changed, and is thus indetermined to the user.
00313 
00314   template <class T, class Signature>
00315     class matrix {
00316 
00317     public:
00318         typedef T value_type;
00319         typedef T scalar_type;
00320         typedef fortran_compatible_tag compatibility;
00321         typedef MMT_compatible_tag MMT_compatibility;
00322         typedef Signature template_signature;
00323 
00324         typedef ptrdiff_t difference_type;
00325         typedef T& reference;
00326         typedef const T& const_reference;
00327         typedef T* storage_iterator;
00328         typedef const T* const_storage_iterator;
00329 
00330         typedef rank2_iterator<T> row_iterator;
00331         typedef rank2_iterator<const T> const_row_iterator;
00332         typedef rank2_iterator<T> column_iterator;
00333         typedef rank2_iterator<const T> const_column_iterator;
00334 
00335         typedef matrix_view<T> matrix_view_type;
00336 
00337     private:
00338         T* array;
00339         int n, m;
00340 
00341     public:
00342         matrix() : array(0) { n = 0; m = 0; }
00343         matrix(int n_, int m_) : array(new T[n_*m_]), n(n_), m(m_) {}
00344         matrix(const matrix& M) : array(new T[M.n*M.m]), n(M.n), m(M.m) {
00345             std::copy(M.array, M.array+n*m, array);
00346         }
00347         ~matrix() { delete[] array; }
00348 
00349         void resize(int i, int j) {
00350             delete[] array;
00351             array = new T[i*j];
00352             n = i; m = j;
00353         }
00354 
00355         int size(int i) const { return i==0? n : m; }
00356         int size0() const { return n; }
00357         int size1() const { return m; }
00358         int size() const { return n*m; }
00359         const T* data() const { return array; }
00360         const T* const_data() const { return array; }
00361         T* data() { return array; }
00362         void modify() { }
00363 
00364         T const& operator()(int i, int j) const
00365             { return array[i+n*j]; }
00366         T& operator()(int i, int j)
00367             { return array[i+n*j]; }
00368 
00369         //  Iterator bounds
00370         storage_iterator storage_begin() { return array; }
00371         storage_iterator storage_end() { return array + n*m; }
00372         const_storage_iterator storage_begin() const { return array; }
00373         const_storage_iterator storage_end() const { return array + n*m; }
00374         int storage_size() const { return n*m; }
00375         row_iterator row_begin() { return row_iterator(array, 1, n, n*m); }
00376         row_iterator row_end() { return row_iterator(array+n, 1, n, n*m); }
00377         column_iterator column_begin()
00378             { return column_iterator(array, n, 1, n); }
00379         column_iterator column_end()
00380             { return column_iterator(array+n*m, n, 1, n); }
00381 
00382         const_row_iterator row_begin() const
00383             { return const_row_iterator(array, 1, n, n*m); }
00384         const_row_iterator row_end() const
00385             { return const_row_iterator(array+n, 1, n, n*m); }
00386         const_column_iterator column_begin() const
00387             { return const_column_iterator(array, n, 1, n); }
00388         const_column_iterator column_end() const
00389             { return const_column_iterator(array+n*m, n, 1, n); }
00390 
00391         row_iterator row_rbegin()
00392             { return row_iterator(array+n-1, -1, n, n*m); }
00393         row_iterator row_rend()
00394             { return row_iterator(array-1, -1, n, n*m); }
00395     };
00396 
00397   //  Note that indexing iterators referring to an non-const hermitian
00398   //  packed matrix violates the iterator requirement that
00399   //  dereferencing it should return a reference.  Instead a proxy is
00400   //  returned in this case.  (With some STL magic, we might define
00401   //  `reference' to this proxy.)
00402 
00403   template <class T>
00404     class matrix<T, hermitian_s> {
00405 
00406       public:
00407         typedef matrix self;
00408         typedef T value_type;
00409         typedef T scalar_type;
00410         typedef fortran_compatible_tag compatibility;
00411         typedef MMT_compatible_tag MMT_compatibility;
00412         typedef hermitian_s template_signature;
00413 
00414         // typedef instance<T>::allocator_type::difference_type
00415         // difference_type; etc., but that is not implemented in class
00416         // instance.
00417         typedef ptrdiff_t difference_type;
00418         typedef T& reference;
00419         typedef const T& const_reference;
00420         typedef T* storage_iterator;
00421         typedef const T* const_storage_iterator;
00422 
00423         typedef indexing_rank2_iterator<self> row_iterator;
00424         typedef indexing_rank2_iterator<self> column_iterator;
00425         typedef indexing_rank2_iterator<const self> const_row_iterator;
00426         typedef indexing_rank2_iterator<const self> const_column_iterator;
00427 
00428       private:
00429         T *array;
00430         int n;
00431 
00432       public:
00433         matrix() : array(0), n(0) {}
00434         matrix(int n_, int m_)
00435             : array(new T[n_*(n_+1)/2]), n(n_) { assert(n_ == m_); }
00436         matrix(const matrix& M) : array(new T[M.size()]), n(M.n) {
00437             std::copy(M.data(), M.data() + M.size(), data());
00438         }
00439         ~matrix() { delete[] array; }
00440 
00441         void resize(int n_) {
00442             delete[] array;
00443             array = new T[n_*(n_+1)/2];
00444             n = n_;
00445         }
00446 
00447         int size(int i) const { return n; }
00448         int size0() const { return n; }
00449         int size1() const { return n; }
00450         int size() const { return n*(n+1)/2; }
00451 
00452         const T* data() const { return array; }
00453         const T* const_data() { return array; }
00454         T* data() { return array; }
00455         void modify() { }
00456         T operator()(int i, int j) const {
00457             if(i <= j) return array[i+j*(j+1)/2];
00458             else       return adj(array[j+i*(i+1)/2]);
00459         }
00460 
00461         class adj_proxy {
00462             T& dest; bool c;
00463         public:
00464             adj_proxy(T& dest_, bool c_) : dest(dest_), c(c_) {}
00465             operator T()
00466                 { if(c) return adj(dest); else return dest; }
00467             adj_proxy operator=(T x)
00468                 { if(c) dest = adj(x); else dest = x; return *this; }
00469             adj_proxy operator+=(T x)
00470                 { if(c) dest += adj(x); else dest += x; return *this; }
00471             adj_proxy operator-=(T x)
00472                 { if(c) dest -= adj(x); else dest -= x; return *this; }
00473             adj_proxy operator*=(T x)
00474                 { if(c) dest *= adj(x); else dest *= x; return *this; }
00475             adj_proxy operator/=(T x)
00476                 { if(c) dest /= adj(x); else dest /= x; return *this; }
00477         };
00478         typedef adj_proxy value_proxy;
00479 
00480         adj_proxy operator()(int i, int j) {
00481             if(i <= j)
00482                 return adj_proxy(array[i+j*(j+1)/2], false);
00483             else
00484                 return adj_proxy(array[j+i*(i+1)/2], true);
00485         }
00486 
00487         const T& upper_triangle(int i, int j) const
00488             { assert(i<=j); return array[i+j*(j+1)/2]; }
00489         T& upper_triangle(int i, int j)
00490             { assert(i<=j); return array[i+j*(j+1)/2]; }
00491 
00492         storage_iterator storage_begin() { return array; }
00493         storage_iterator storage_end() { return array + n*(n+1)/2; }
00494         const_storage_iterator storage_begin() const { return array; }
00495         const_storage_iterator storage_end() const
00496             { return array + n*(n+1)/2; }
00497         int storage_size() const { return n*(n+1)/2; }
00498 
00499         row_iterator row_begin() { return row_iterator(*this, n, 0); }
00500         row_iterator row_end() { return row_iterator(*this, n, n); }
00501         column_iterator column_begin() { return column_iterator(*this,n,0); }
00502         column_iterator column_end() { return column_iterator(*this,n,n); }
00503         const_row_iterator row_begin() const
00504             { return const_row_iterator(*this, n, 0); }
00505         const_row_iterator row_end() const
00506             { return const_row_iterator(*this, n, n); }
00507         const_column_iterator column_begin() const
00508             { return const_column_iterator(*this,n,0); }
00509         const_column_iterator column_end() const
00510             { return const_column_iterator(*this,n,n); }
00511     };
00512 
00513   template <class T>
00514     class mvector {
00515     public:
00516         typedef T value_type;
00517         typedef T scalar_type;
00518         typedef T* iterator;
00519         typedef const T* const_iterator;
00520         typedef iterator column_iterator;
00521         typedef const_iterator const_column_iterator;
00522         typedef iterator row_iterator;
00523         typedef const_iterator const_row_iterator;
00524         typedef iterator diagonal_iterator;
00525         typedef const_iterator const_diagonal_iterator;
00526         typedef int size_type;
00527 
00528     private:
00529         size_type n;
00530         T* array;
00531 
00532     public:
00533         mvector() : n(0), array(new T[1]) {}
00534         explicit mvector(int n_) : n(n_), array(new T[n_+!n_]) {}
00535         mvector(int n_, int m_) : n(n_), array(new T[n_+!n_])
00536             { assert(m_==1); }
00537         mvector(const mvector& x) : n(x.n), array(new T[x.n+!x.n])
00538             { std::copy(x.array, x.array + n, array); }
00539         ~mvector() { delete[] array; }
00540         mvector& operator=(mvector const& rhs) {
00541             delete[] array;
00542             n = rhs.n;
00543             array = new T[rhs.n+!rhs.n];
00544             std::copy(rhs.array, rhs.array + n, array);
00545             return *this;
00546         }
00547 
00548         void resize(int i) {
00549             delete[] array;
00550             array = new T[i];
00551             n = i;
00552         }
00553         bool empty() const { return n == 0; }
00554 
00555         int size() const { return n; }
00556         int size(int i) const { return i==0? n : 1; }
00557         int size0() const { return n; }
00558         int size1() const { return 1; }
00559         const value_type* const_data() const { return array; }
00560         const value_type* data() const { return array; }
00561         value_type* data() { return array; }
00562         void modify() { }
00563 
00564         T const& operator()(int i, int j) const
00565             { assert(j==0); assert(i >= 0 && i < n); return array[i]; }
00566         T& operator()(int i, int j)
00567             { assert(j==0); assert(i >= 0 && i < n); return array[i]; }
00568         T const& operator()(int i) const
00569             { assert(i>=0 && i<n); return array[i]; }
00570         T& operator()(int i) { assert(i>=0 && i<n); return array[i]; }
00571         T const& operator[](int i) const
00572             { assert(i>=0 && i<n); return array[i]; }
00573         T& operator[](int i) { assert(i>=0 && i<n); return array[i]; }
00574 
00575         const_iterator begin() const { return array; }
00576         const_iterator end() const { return array + n; }
00577         iterator begin() { return array; }
00578         iterator end() { return array + n; }
00579 
00580         const_column_iterator column_begin(int i = 0) const
00581             { return begin(); }
00582         const_column_iterator column_end(int i = 0) const
00583             { return end(); }
00584         column_iterator column_begin(int i = 0) { return begin(); }
00585         column_iterator column_end(int i = 0) { return end(); }
00586 
00587         const_row_iterator row_begin(int i) const
00588             { return array + i; }
00589         const_row_iterator row_end(int i) const
00590             { return array + i + 1; }
00591         row_iterator row_begin(int i)
00592             { return array + i; }
00593         row_iterator row_end(int i)
00594             { return array + i + 1; }
00595 
00596         const_diagonal_iterator diagonal_begin(int i) const
00597             { return row_begin(-i); }
00598         const_diagonal_iterator diagonal_end(int i) const
00599             { return row_end(-i); }
00600         diagonal_iterator diagonal_begin(int i) { return row_begin(-i); }
00601         diagonal_iterator diagonal_end(int i) { return row_end(-i); }
00602 
00603         mvector& operator+=(mvector const& rhs) {
00604             assert(n == rhs.n);
00605             for(int i = 0; i < n; ++i) array[i] += rhs.array[i];
00606             return *this;
00607         }
00608         mvector& operator-=(mvector const& rhs) {
00609             assert(n == rhs.n);
00610             for(int i = 0; i < n; ++i) array[i] -= rhs.array[i];
00611             return *this;
00612         }
00613         template<typename U>
00614         mvector& operator*=(U const& rhs) {
00615             for(int i = 0; i < n; ++i) array[i] *= rhs;
00616             return *this;
00617         }
00618         template<typename U>
00619         mvector& operator/=(U const& rhs) {
00620             for(int i = 0; i < n; ++i) array[i] /= rhs;
00621             return *this;
00622         }
00623     };
00624 }
00625 // namespace math {
00626 //   template<typename T>
00627 //     struct scalar_type_< mvector<T> > { typedef T eval; };
00628 //   template<typename T, typename Sig>
00629 //     struct scalar_type_< matrix<T, Sig> > { typedef T eval; };
00630 // }
00631 }
00632 ///\endif
00633 
00634 #endif