nemo/Matrix.h
Go to the documentation of this file.00001 #ifndef _NEMO_MATRIX_H_ 00002 #define _NEMO_MATRIX_H_ 00003 00004 00005 #include "IntrusiveCowSupport.h" 00006 #include "IntrusiveCounted.h" 00007 #include "Vector.h" 00008 #include "Collection.h" 00009 #include "detail/nemo-cpp11.h" 00010 #include "detail/nemo-errors.h" 00011 00012 #include <Eigen/Core> 00013 00014 #include <stdexcept> 00015 #include <sstream> 00016 00017 00018 namespace nemo { 00019 00020 00021 class MatrixDimension; 00022 static MatrixDimension dim(unsigned int rows, unsigned int cols); 00023 00028 class MatrixDimension { 00029 public: 00033 inline const unsigned int& rows() const {return numRows;} 00037 inline const unsigned int& cols() const {return numCols;} 00038 00039 friend MatrixDimension dim(unsigned int rows, unsigned int cols); 00040 MatrixDimension(const MatrixDimension &o) : numRows(o.numRows), numCols(o.numCols) {} 00041 MatrixDimension& operator=(const MatrixDimension& o){numRows=o.numRows; numCols=o.numCols; return *this;} 00042 bool operator==(const MatrixDimension &o)const{return numCols==o.numCols&&numRows==o.numRows;} 00043 bool operator!=(const MatrixDimension &o)const{return !(*this==o);} 00044 friend std::ostream& operator<<(std::ostream &o, const MatrixDimension &d){ 00045 o<<std::string("dim(")<<(d.numRows)<<std::string(",")<<(d.numRows)<<std::string(")"); 00046 return o; 00047 } 00048 private: 00049 MatrixDimension(unsigned int r, unsigned int c) : numRows(r), numCols(c) {} 00050 unsigned int numRows; 00051 unsigned int numCols; 00052 }; 00053 00057 static inline MatrixDimension dim(unsigned int rows, unsigned int cols){ 00058 return MatrixDimension(rows, cols); 00059 } 00060 00061 00062 // forward declaration for friends in MatrixMath.h 00063 template <class T> 00064 class Matrix; 00065 #ifndef NEMO_PARSED_BY_DOXYGEN 00066 template <typename T> 00067 Matrix<T> matMatT(const Matrix<T>&); 00068 template <typename T> 00069 void _nemo_impl_matMatTRecursive(const Matrix<T>&,Matrix<T> &,unsigned int, 00070 unsigned int,unsigned int, unsigned int,unsigned int,bool); 00071 template <typename T> 00072 void _nemo_impl_matMatTBlockRecursive(const Matrix<T>&,Matrix<T> &,unsigned int, 00073 unsigned int,unsigned int,unsigned int,unsigned int,unsigned int,bool); 00074 template <typename T> 00075 void _nemo_impl_matMatTBlockIterative(const Matrix<T> &m, Matrix<T> &result); 00076 #endif 00077 00078 00079 #ifndef NEMO_PARSED_BY_DOXYGEN 00080 #ifdef NEMO_ABORT_ON_ERROR 00081 #define _NEMO_MATRIX_CHECK_INDEX(rowIndex,colIndex) \ 00082 if(rows() <= rowIndex || cols() <= colIndex){ \ 00083 abort(); \ 00084 } 00085 #else 00086 #define _NEMO_MATRIX_CHECK_INDEX(rowIndex,colIndex) \ 00087 if(rows() <= rowIndex){ \ 00088 throw std::invalid_argument("Matrix: row-index out of range"); \ 00089 } \ 00090 if(cols() <= colIndex){ \ 00091 throw std::invalid_argument("Matrix: column-index out of range"); \ 00092 } 00093 #endif 00094 #endif 00095 00116 template <class T> 00117 class Matrix : public Collection<T> { 00118 private: 00119 #ifndef NEMO_PARSED_BY_DOXYGEN 00120 //friend class MatrixMath<T>; 00121 //friend Matrix<T> pseudoInverse(Matrix<T>&); 00122 00123 //template <typename D> 00124 //friend Matrix<D> matMatT(const Matrix<D>&); 00125 //template <typename D> 00126 //friend void _nemo_impl_matMatTRecursive(const Matrix<D>&,Matrix<D> &,unsigned int, 00127 // unsigned int,unsigned int, unsigned int,unsigned int,bool); 00128 friend Matrix<T> matMatT<T>(const Matrix<T>&); 00129 friend void _nemo_impl_matMatTRecursive<T>(const Matrix<T>&,Matrix<T> &,unsigned int, 00130 unsigned int,unsigned int, unsigned int,unsigned int,bool); 00131 friend void _nemo_impl_matMatTBlockRecursive<T>(const Matrix<T>&,Matrix<T> &,unsigned int, 00132 unsigned int,unsigned int, unsigned int,unsigned int,unsigned int,bool); 00133 friend void _nemo_impl_matMatTBlockIterative<T>(const Matrix<T> &m, Matrix<T> &result); 00134 #endif 00135 00136 #ifndef NEMO_PARSED_BY_DOXYGEN 00137 class MatrixData : public IntrusiveCounted<MatrixData > { 00138 public: 00139 MatrixData(const MatrixDimension &dimP) 00140 : IntrusiveCounted<MatrixData >(), dim(dimP), mat() 00141 { 00142 // note eigen does only allow 0x0, but not 0xN matrices ^^ 00143 if(dim.rows()>0&&dim.cols()>0){ 00144 mat = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>((int)dim.rows(),(int)dim.cols()); 00145 } 00146 } 00147 MatrixData(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &matP) 00148 : IntrusiveCounted<MatrixData >(), dim(nemo::dim((unsigned int)matP.rows(), (unsigned int)matP.cols())), mat(matP) 00149 { 00150 } 00151 00152 MatrixData(const MatrixData &md) 00153 : IntrusiveCounted<MatrixData >(), dim(md.dim), mat(md.mat) 00154 {} 00155 00156 MatrixData& operator=(const MatrixData &r){ 00157 if(this != &r){ 00158 dim = r.dim; 00159 mat = r.mat; 00160 } 00161 return *this; 00162 } 00163 00164 ~MatrixData(){ 00165 } 00166 public: 00167 mutable MatrixDimension dim; 00168 mutable Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> mat; 00169 }; 00170 #endif 00171 00172 public: 00173 // ============================================================================= 00174 // Construction 00175 // ============================================================================= 00176 00180 explicit Matrix() 00181 : data(new MatrixData(dim(0,0))) 00182 {} 00183 00191 explicit Matrix(const MatrixDimension& dim) 00192 : data(new MatrixData(dim)) 00193 { 00194 for(unsigned int r=0; r<dim.rows(); r++){ 00195 for(unsigned int c=0; c<dim.cols(); c++){ 00196 data->mat((int)r,(int)c) = T(); 00197 } 00198 } 00199 } 00200 00208 explicit Matrix(const MatrixDimension &dim, const T &val) 00209 : data(new MatrixData(dim)) 00210 { 00211 for(unsigned int r=0; r<dim.rows(); r++){ 00212 for(unsigned int c=0; c<dim.cols(); c++){ 00213 // FIXME eigen is indexed with signed ints 00214 // => check maxDimension to assure no sign-conversion 00215 data->mat((int)r,(int)c) = val; 00216 } 00217 } 00218 } 00219 00230 template <typename generatorType> 00231 explicit Matrix(const MatrixDimension &dim, generatorType gen) 00232 : data(new MatrixData(dim)) 00233 { 00234 for(unsigned int r=0; r<dim.rows(); r++){ 00235 for(unsigned int c=0; c<dim.cols(); c++){ 00236 data->mat((int)r,(int)c) = gen(); 00237 } 00238 } 00239 } 00240 00241 #ifdef _NEMO_HAS_INIT_LISTS_ 00242 Matrix(const std::initializer_list<std::initializer_list<T>> &init) 00243 : data(new MatrixData(dim(0,0))) 00244 { 00245 unsigned int numRows = init.size(); 00246 unsigned int numCols = 0; 00247 for(const std::initializer_list<T> *p = init.begin (); p != init.end (); ++p){ 00248 if( (*p).size() > numCols ){ 00249 numCols = (*p).size(); 00250 } 00251 } 00252 data = IntrusiveCowSupport<MatrixData>(new MatrixData(dim(numRows,numCols))); 00253 if(numRows == 0 || numCols == 0){ 00254 return; // no internal Eigen matrix created 00255 } 00256 unsigned int rowIndex = 0; 00257 for(const std::initializer_list<T> *rowP = init.begin(); rowP != init.end(); ++rowP){ 00258 unsigned int colIndex = 0; 00259 for(const T *val = (*rowP).begin(); val != (*rowP).end(); ++val){ 00260 data->mat((int)rowIndex,(int)(colIndex++)) = *val; 00261 } 00262 while(colIndex <= numCols-1){ 00263 data->mat((int)rowIndex,(int)(colIndex++)) = T(); 00264 } 00265 rowIndex++; 00266 } 00267 } 00268 #endif 00269 00278 static Matrix<T> diagonal(const MatrixDimension &dim, const T &diagValue){ 00279 Matrix<T> matrix(dim); 00280 unsigned int diagCount = (dim.rows()<dim.cols()) ? dim.rows() : dim.cols(); 00281 00282 for(unsigned int i=0; i<diagCount; i++){ 00283 matrix.data->mat((int)i,(int)i) = diagValue; 00284 } 00285 return matrix; 00286 } 00287 00296 static Matrix<T> diagonal(const MatrixDimension &dim, const MathVector<T> &diag){ 00297 unsigned int diagCount = (dim.rows()<dim.cols()) ? dim.rows() : dim.cols(); 00298 if(diagCount != diag.dimension()){ 00299 _NEMO_ARGUMENT_ERROR("dimension mismatch in diagonal matrix contruction"); 00300 } 00301 Matrix<T> matrix(dim); 00302 for(unsigned int i=0; i<diagCount; i++){ 00303 matrix.data->mat((int)i,(int)i) = diag[i]; 00304 } 00305 return matrix; 00306 } 00315 static Matrix<T> diagonal(const MathVector<T> &diag){ 00316 Matrix<T> matrix(dim(diag.dimension(), diag.dimension())); 00317 for(unsigned int i=0; i<diag.dimension(); i++){ 00318 matrix.data->mat((int)i,(int)i) = diag[i]; 00319 } 00320 return matrix; 00321 } 00322 00323 // cppcheck-suppress operatorEq 00324 Matrix<T>& operator=(const Matrix<T> &m){ 00325 this->data = m.data; 00326 return *this; 00327 } 00328 00329 //private: 00330 explicit Matrix(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &m) 00331 : data(new MatrixData(m)) 00332 {} 00333 //explicit Matrix(Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> m) 00334 // : data(new MatrixData(m)) 00335 //{} 00336 // ============================================================================= 00337 // Access 00338 // ============================================================================= 00339 public: 00343 MatrixDimension dimension() const { 00344 return data->dim; 00345 } 00346 00350 unsigned int rows() const { 00351 return data->dim.rows(); 00352 } 00353 00357 unsigned int cols() const { 00358 return data->dim.cols(); 00359 } 00360 00361 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> eigen() const { 00362 return data->mat; 00363 } 00364 00368 const T& operator()(const unsigned int rowIndex,const unsigned int colIndex) const{ 00369 _NEMO_MATRIX_CHECK_INDEX(rowIndex, colIndex); 00370 return data->mat((int)rowIndex,(int)colIndex); 00371 } 00372 00373 #ifndef NEMO_PARSED_BY_DOXYGEN 00374 class ReadWriteRef { 00375 friend class Matrix; 00376 private: 00377 ReadWriteRef(IntrusiveCowSupport<MatrixData > &d, const unsigned int r, const unsigned int c) 00378 : data(d), rowIndex(r), colIndex(c) {} 00379 public: 00380 ReadWriteRef& operator=(const T &value){ 00381 // write access 00382 data.detach(); 00383 data->mat((int)rowIndex,(int)colIndex) = value; 00384 return *this; 00385 } 00386 ReadWriteRef& operator=(ReadWriteRef &value){ 00387 // write access 00388 data.detach(); 00389 data->mat((int)rowIndex,(int)colIndex) = (T)value; 00390 return *this; 00391 } 00392 00393 operator const T() const{return data->mat((int)rowIndex,(int)colIndex);} 00394 00395 #define _nemo_matrix_ref_define_compound_(OP) \ 00396 ReadWriteRef& operator OP (const T& value){\ 00397 data.detach();\ 00398 data->mat((int)rowIndex,(int)colIndex) OP value;\ 00399 return *this;\ 00400 } 00401 _nemo_matrix_ref_define_compound_(+=) 00402 _nemo_matrix_ref_define_compound_(-=) 00403 _nemo_matrix_ref_define_compound_(*=) 00404 _nemo_matrix_ref_define_compound_(/=) 00405 _nemo_matrix_ref_define_compound_(%=) 00406 _nemo_matrix_ref_define_compound_(&=) 00407 _nemo_matrix_ref_define_compound_(|=) 00408 _nemo_matrix_ref_define_compound_(^=) 00409 _nemo_matrix_ref_define_compound_(<<=) 00410 _nemo_matrix_ref_define_compound_(>>=) 00411 00412 #define _nemo_matrix_ref_define_compare_(OP) \ 00413 bool operator OP (const typename Matrix<T>::ReadWriteRef& rhs) const{\ 00414 return data->mat((int)rowIndex,(int)colIndex) OP rhs.data->mat((int)rhs.rowIndex,(int)rhs.colIndex);\ 00415 } 00416 _nemo_matrix_ref_define_compare_(==) 00417 _nemo_matrix_ref_define_compare_(!=) 00418 _nemo_matrix_ref_define_compare_(<=) 00419 _nemo_matrix_ref_define_compare_(<) 00420 _nemo_matrix_ref_define_compare_(>=) 00421 _nemo_matrix_ref_define_compare_(>) 00422 00423 // this is a 'specialization' of std::swap that can be found by argument-dependent lookup. 00424 // this one is necessary for STL write operations to work correctly, since the standard swap 00425 // implementation does not work correctly on ReadWriteRef, since it already has reference semantics 00426 friend inline void swap( typename nemo::Matrix<T>::ReadWriteRef a, typename nemo::Matrix<T>::ReadWriteRef b ){ 00427 a.data.detach(); 00428 b.data.detach(); 00429 using std::swap; 00430 swap( 00431 a.data->mat((int)a.rowIndex,(int)a.colIndex), 00432 b.data->mat((int)b.rowIndex,(int)b.colIndex) 00433 ); 00434 } 00435 friend inline T max( const typename nemo::Matrix<T>::ReadWriteRef& a, const typename nemo::Matrix<T>::ReadWriteRef& b ){ 00436 using std::max; 00437 return max( 00438 a.data->mat((int)a.rowIndex,(int)a.colIndex), 00439 b.data->mat((int)b.rowIndex,(int)b.colIndex) 00440 ); 00441 } 00442 friend inline T max( const typename nemo::Matrix<T>::ReadWriteRef& a, T b ){ 00443 using std::max; 00444 return max( a.data->mat((int)a.rowIndex,(int)a.colIndex), b ); 00445 } 00446 friend inline T max( T a, const typename nemo::Matrix<T>::ReadWriteRef& b ){ 00447 using std::max; 00448 return max( a, b.data->mat((int)b.rowIndex,(int)b.colIndex) ); 00449 } 00450 friend inline T min( const typename nemo::Matrix<T>::ReadWriteRef& a, const typename nemo::Matrix<T>::ReadWriteRef& b ){ 00451 using std::min; 00452 return min( 00453 a.data->mat((int)a.rowIndex,(int)a.colIndex), 00454 b.data->mat((int)b.rowIndex,(int)b.colIndex) 00455 ); 00456 } 00457 friend inline T min( const typename nemo::Matrix<T>::ReadWriteRef& a, T b ){ 00458 using std::min; 00459 return min( a.data->mat((int)a.rowIndex,(int)a.colIndex), b ); 00460 } 00461 friend inline T min( T a, const typename nemo::Matrix<T>::ReadWriteRef& b ){ 00462 using std::min; 00463 return min( a, b.data->mat((int)b.rowIndex,(int)b.colIndex) ); 00464 } 00465 00466 private: 00467 IntrusiveCowSupport<MatrixData > &data; 00468 unsigned int rowIndex; 00469 unsigned int colIndex; 00470 }; 00471 #endif 00472 00477 ReadWriteRef operator()(const unsigned int rowIndex, const unsigned int colIndex){ 00478 _NEMO_MATRIX_CHECK_INDEX(rowIndex,colIndex); 00479 return ReadWriteRef(data, rowIndex, colIndex); 00480 } 00481 00482 00486 Matrix<T> subMatrix(const unsigned int fromRow, const unsigned int fromCol, const MatrixDimension &d) const { 00487 checkSubsumes(fromRow, fromCol, d); 00488 00489 // prevent Eigen from going mad about 0-dimensional matrices 00490 if( d.cols()==0 || d.rows()==0 ){ 00491 return Matrix<T>(d); 00492 } 00493 00494 Matrix<T> sub(d); 00495 sub.data->mat = this->data->mat.block((int)fromRow,(int)fromCol,(int)d.rows(),(int)d.cols()); 00496 return sub; 00497 } 00498 00503 void setSubMatrix(const unsigned int fromRow, const unsigned int fromCol, const Matrix<T> &m){ 00504 MatrixDimension d = m.dimension(); 00505 checkSubsumes(fromRow, fromCol, d); 00506 00507 // prevent Eigen from going mad about 0-dimensional matrices 00508 if( d.cols()==0 || d.rows()==0 ){ 00509 return; 00510 } 00511 00512 this->data->mat.block((int)fromRow,(int)fromCol,(int)d.rows(),(int)d.cols()) 00513 = m.data->mat; 00514 } 00515 00517 MathVector<T> col(unsigned int colIndex){ 00518 if(cols() <= colIndex){ 00519 _NEMO_ARGUMENT_ERROR("Matrix::col(uint): column index out of range"); 00520 } 00521 MathVector<T> c(dim(this->rows())); 00522 for(unsigned int i=0; i<rows(); i++){ 00523 c[i] = data->mat((int)i, (int)colIndex); 00524 } 00525 return c; 00526 } 00527 00529 void setCol(unsigned int colIndex, const MathVector<T> &vec){ 00530 if(cols() <= colIndex){ 00531 _NEMO_ARGUMENT_ERROR("Matrix::setCol(uint,MathVector): column index out of range"); 00532 } 00533 if(vec.dimension() != this->rows()){ 00534 _NEMO_ARGUMENT_ERROR("Matrix::setCol(uint,MathVector): dimension mismatch"); 00535 } 00536 00537 for(unsigned int i=0; i<rows(); i++){ 00538 data->mat((int)i, (int)colIndex) = vec[i]; 00539 } 00540 } 00541 00543 VectorTranspose<T> row(unsigned int rowIndex){ 00544 if(rows() <= rowIndex){ 00545 _NEMO_ARGUMENT_ERROR("Matrix::row(uint): row index out of range"); 00546 } 00547 MathVector<T> r(dim(this->cols())); 00548 for(unsigned int i=0; i<cols(); i++){ 00549 r[i] = data->mat((int)rowIndex, (int)i); 00550 } 00551 return r.transpose(); 00552 } 00553 00555 void setRow(unsigned int rowIndex, const VectorTranspose<T> &vec){ 00556 if(rows() <= rowIndex){ 00557 _NEMO_ARGUMENT_ERROR("Matrix::setRow(uint,MathVector): row index out of range"); 00558 } 00559 if(vec.dimension() != this->cols()){ 00560 _NEMO_ARGUMENT_ERROR("Matrix::setRow(uint,MathVector): dimension mismatch"); 00561 } 00562 00563 for(unsigned int i=0; i<cols(); i++){ 00564 data->mat((int)rowIndex, (int)i) = vec[i]; 00565 } 00566 } 00567 00568 // ============================================================================= 00569 // Comparison 00570 // ============================================================================= 00575 bool operator==(const Matrix<T> &o) const{ 00576 if(data->dim != o.data->dim){ 00577 return false; 00578 } 00579 for(unsigned int r=0; r<rows(); r++){ 00580 for(unsigned int c=0; c<cols(); c++){ 00581 if(data->mat((int)r,(int)c) != o.data->mat((int)r,(int)c)){ 00582 return false; 00583 } 00584 } 00585 } 00586 // else 00587 return true; 00588 } 00589 00593 bool operator!=(const Matrix<T> &r) const{ 00594 return !operator==(r); 00595 } 00596 00597 // ============================================================================= 00598 // Collection API 00599 // ============================================================================= 00600 00601 unsigned int numElements() const { 00602 return rows() * cols(); 00603 } 00604 T getElement(unsigned int elementIndex) const{ 00605 if(numElements() <= elementIndex){ 00606 _NEMO_ARGUMENT_ERROR("element index out of range"); 00607 } 00608 const unsigned int rowIndex = elementIndex / cols(); 00609 return (*this)(rowIndex,elementIndex-rowIndex*cols()); 00610 } 00611 void setElement(unsigned int elementIndex, T value){ 00612 if(numElements() <= elementIndex){ 00613 _NEMO_ARGUMENT_ERROR("element index out of range"); 00614 } 00615 const unsigned int rowIndex = elementIndex / cols(); 00616 (*this)(rowIndex,elementIndex-rowIndex*cols()) = value; 00617 } 00618 00619 void transform(Mapping<T,T> m){ 00620 for(unsigned int i=0; i<rows(); i++){ 00621 for(unsigned int k=0; k<cols(); k++){ 00622 (*this)(i,k) = m((*this)(i,k)); 00623 } 00624 } 00625 } 00626 00627 template <typename OutType> 00628 Matrix<OutType> map(Mapping<T,OutType> m){ 00629 Matrix<OutType> out(this->dimension()); 00630 for(unsigned int i=0; i<rows(); i++){ 00631 for(unsigned int k=0; k<cols(); k++){ 00632 out(i,k) = m((*this)(i,k)); 00633 } 00634 } 00635 return out; 00636 } 00637 00638 // ============================================================================= 00639 // STL compliant iterator interface 00640 // ============================================================================= 00641 00642 #ifndef NEMO_PARSED_BY_DOXYGEN 00643 class Iterator{ 00644 private: 00645 Iterator(const Matrix<T> *matP, int indexP) : mat(matP), index(indexP), rwRef(const_cast<Matrix<T>*>(matP)->data,0,0) {} 00646 public: 00647 Iterator(const Iterator& it) 00648 : mat(it.mat), index(it.index), 00649 rwRef(const_cast<Matrix<T>*>(it.mat)->data,(unsigned int)it.rwRef.rowIndex,(unsigned int)it.rwRef.colIndex) 00650 {} 00651 Iterator& operator=(const Iterator &it){ 00652 if(this==&it){return *this;} 00653 this->mat = it.mat; 00654 this->index = it.index; 00655 this->rwRef.data = it.rwRef.data; 00656 this->rwRef.rowIndex = it.rwRef.rowIndex; 00657 this->rwRef.colIndex = it.rwRef.colIndex; 00658 return *this; 00659 } 00660 00661 typedef int difference_type; 00662 typedef T value_type; 00663 typedef ReadWriteRef* pointer; 00664 typedef ReadWriteRef reference; 00665 typedef std::random_access_iterator_tag iterator_category; 00666 00667 ReadWriteRef& operator*(){ 00668 if(index<0 || (int)mat->rows()*(int)mat->cols()<=index){_NEMO_LOGIC_ERROR("Iterator out of range");} 00669 rwRef.colIndex = (unsigned int)(index%(int)mat->cols()); 00670 rwRef.rowIndex = (unsigned int)(index/(int)mat->cols()); 00671 return rwRef; 00672 //return (*const_cast<MathVector<T>*>(vec))[(unsigned int)index]; 00673 } 00674 T operator*() const{ 00675 if(index<0 || (int)mat->rows()*(int)mat->cols()<=index){_NEMO_LOGIC_ERROR("Iterator out of range");} 00676 return mat->getElement((unsigned int)index); 00677 } 00678 ReadWriteRef& operator[](difference_type offset){ 00679 if(index+offset<0 || (int)mat->rows()*(int)mat->cols()<=index+offset){_NEMO_LOGIC_ERROR("Iterator out of range");} 00680 rwRef.colIndex = (unsigned int)(index%(int)mat->cols()); 00681 rwRef.rowIndex = (unsigned int)(index/(int)mat->cols()); 00682 rwRef.index = (unsigned int)index+offset; 00683 return rwRef; 00684 //return (*const_cast<MathVector<T>*>(vec))[(unsigned int)index+offset]; 00685 } 00686 T operator[](difference_type offset) const{ 00687 if(index+offset<0 || (int)mat->rows()*(int)mat->cols()<=index+offset){_NEMO_LOGIC_ERROR("Iterator out of range");} 00688 return mat->getElement((unsigned int)index+offset); 00689 } 00690 Iterator& operator++() const{ 00691 index++; 00692 return const_cast<Iterator&>(*this); 00693 } 00694 Iterator operator++(int) const{ 00695 Iterator copy = *this; 00696 index++; 00697 return copy; 00698 } 00699 Iterator& operator--() const{ 00700 index--; 00701 return const_cast<Iterator&>(*this); 00702 } 00703 Iterator operator--(int) const{ 00704 Iterator copy = *this; 00705 index--; 00706 return copy; 00707 } 00708 Iterator operator+(difference_type off) const { 00709 Iterator it = *this; 00710 it.index += off; 00711 return it; 00712 } 00713 Iterator operator-(difference_type off) const { 00714 Iterator it = *this; 00715 it.index -= off; 00716 return it; 00717 } 00718 difference_type operator-(const Iterator& rhs) const { 00719 return index - rhs.index; 00720 } 00721 bool operator==(const Iterator &it) const{ 00722 return (mat==it.mat) && (index==it.index); 00723 } 00724 bool operator!=(const Iterator &it) const{ 00725 return !(*this==it); 00726 } 00727 bool operator<=(const Iterator &it) const{ 00728 return index<=it.index; 00729 } 00730 bool operator<(const Iterator &it) const{ 00731 return index<it.index; 00732 } 00733 bool operator>=(const Iterator &it) const{ 00734 return index>=it.index; 00735 } 00736 bool operator>(const Iterator &it) const{ 00737 return index>it.index; 00738 } 00739 private: 00740 friend class Matrix<T>; 00741 const Matrix<T> *mat; 00742 mutable int index; 00743 mutable ReadWriteRef rwRef; 00744 }; 00745 #endif 00746 00750 Iterator begin(){ 00751 return Iterator(this,0); 00752 } 00756 Iterator end(){ 00757 return Iterator(this,(int)rows()*(int)cols()); 00758 } 00759 00760 const Iterator begin() const{ 00761 return Iterator(this,0); 00762 } 00763 const Iterator end() const{ 00764 return Iterator(this,(int)rows()*(int)cols()); 00765 } 00766 00767 00768 // ============================================================================= 00769 // Unary Matrix operations 00770 // ============================================================================= 00771 00772 Matrix<T> transpose() const { 00773 Matrix<T> t(dim(this->cols(), this->rows())); 00774 t.data->mat = this->data->mat.transpose(); 00775 return t; 00776 } 00777 00778 00779 // ============================================================================= 00780 // Matrix-Matrix operations 00781 // ============================================================================= 00782 00786 const Matrix<T> attachRows(const Matrix<T> &o) const { 00787 checkColumnDimensionMatch(o); 00788 00789 Matrix<T> newMat(dim(rows()+o.rows(), cols())); 00790 if(rows()>0 && cols()>0){ 00791 newMat.data->mat.block(0,0,(int)rows(),(int)cols()) = this->data->mat; 00792 } 00793 if(o.rows()>0 && o.cols()>0){ 00794 newMat.data->mat.block((int)rows(),0,(int)o.rows(),(int)o.cols()) = o.data->mat; 00795 } 00796 return newMat; 00797 } 00798 00802 const Matrix<T> attachCols(const Matrix<T> &o) const { 00803 checkRowDimensionMatch(o); 00804 00805 Matrix<T> newMat(dim(rows(), cols()+o.cols())); 00806 if(rows()>0 && cols()>0){ 00807 newMat.data->mat.block(0,0,(int)rows(),(int)cols()) = this->data->mat; 00808 } 00809 if(o.rows()>0 && o.cols()>0){ 00810 newMat.data->mat.block(0,(int)cols(),(int)o.rows(),(int)o.cols()) = o.data->mat; 00811 } 00812 return newMat; 00813 } 00814 00815 private: 00816 #ifndef NEMO_PARSED_BY_DOXYGEN 00817 class MatrixMult{ 00818 public: 00819 static void evaluate(const Matrix<T> &m1, const Matrix<T> &m2, Matrix<T> &result){ 00820 if(m1.rows()!=0 && m2.cols()!=0 && m1.cols()==0){ 00821 result.data->mat = Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Zero(m1.rows(), m2.cols()); 00822 } 00823 else{ 00824 result.data->mat = m1.data->mat * m2.data->mat; 00825 } 00826 } 00827 }; 00828 class MatrixAdd{ 00829 public: 00830 static void evaluate(const Matrix<T> &m1, const Matrix<T> &m2, Matrix<T> &result){ 00831 result.data->mat = m1.data->mat + m2.data->mat; 00832 } 00833 }; 00834 class MatrixDiff{ 00835 public: 00836 static void evaluate(const Matrix<T> &m1, const Matrix<T> &m2, Matrix<T> &result){ 00837 result.data->mat = m1.data->mat - m2.data->mat; 00838 } 00839 }; 00840 class MatrixCWiseMult{ 00841 public: 00842 static void evaluate(const Matrix<T> &m1, const Matrix<T> &m2, Matrix<T> &result){ 00843 result.data->mat = m1.data->mat.cwise() * m2.data->mat; 00844 } 00845 }; 00846 class MatrixCWiseDiv{ 00847 public: 00848 static void evaluate(const Matrix<T> &m1, const Matrix<T> &m2, Matrix<T> &result){ 00849 result.data->mat = m1.data->mat.cwise() / m2.data->mat; 00850 } 00851 }; 00852 class MatrixCWiseMax{ 00853 public: 00854 static void evaluate(const Matrix<T> &m1, const Matrix<T> &m2, Matrix<T> &result){ 00855 result.data->mat = m1.data->mat.cwise().max( m2.data->mat ); 00856 } 00857 }; 00858 class MatrixCWiseMin{ 00859 public: 00860 static void evaluate(const Matrix<T> &m1, const Matrix<T> &m2, Matrix<T> &result){ 00861 result.data->mat = m1.data->mat.cwise().min( m2.data->mat ); 00862 } 00863 }; 00864 #endif 00865 00866 public: 00870 Matrix<T> operator*(const Matrix<T> &r) const{ 00871 if(this->cols() != r.rows()){ 00872 _NEMO_ARGUMENT_ERROR("Matrix multiplication: dimension mismatch"); 00873 } 00874 // FIXME [n x 0]*[0 x m] must not evaluate to Eigen::mat [0 x 0] 00875 Matrix<T> result(dim(this->rows(), r.cols())); 00876 MatrixMult::evaluate(*this, r, result); 00877 return result; 00878 } 00883 Matrix<T>& operator*=(const Matrix<T> &r){ 00884 if(this->cols() != r.rows()){ 00885 _NEMO_ARGUMENT_ERROR("Matrix multiplication: dimension mismatch"); 00886 } 00887 data.detach(); 00888 data->dim = dim(this->rows(), r.cols()); 00889 //MatrixDimension newDim = dim(this->rows(), r.cols()); 00890 //data->dim = newDim; 00891 MatrixMult::evaluate(*this, r, *this); 00892 return *this; 00893 } 00897 Matrix<T> operator+(const Matrix<T> &r) const{ 00898 checkDimensionMatch(r); 00899 Matrix<T> result(dimension()); 00900 MatrixAdd::evaluate(*this, r, result); 00901 return result; 00902 } 00906 Matrix<T>& operator+=(const Matrix<T> &r){ 00907 checkDimensionMatch(r); 00908 data.detach(); 00909 MatrixAdd::evaluate(*this, r, *this); 00910 return *this; 00911 } 00915 Matrix<T> operator-(const Matrix<T> &r) const{ 00916 checkDimensionMatch(r); 00917 Matrix<T> result(dimension()); 00918 MatrixDiff::evaluate(*this, r, result); 00919 return result; 00920 } 00924 Matrix<T>& operator-=(const Matrix<T> &r){ 00925 checkDimensionMatch(r); 00926 data.detach(); 00927 MatrixDiff::evaluate(*this, r, *this); 00928 return *this; 00929 } 00930 00931 // ============================================================================= 00932 // Matrix-Matrix component-wise operations 00933 // ============================================================================= 00934 00939 class CWise { 00940 friend class Matrix; 00941 private: 00942 CWise(Matrix<T> &m) : mat(m) {} 00943 CWise(const CWise &o) : mat(o.mat) {} 00944 CWise& operator=(const CWise &o){this->mat=o.mat; return *this;} 00945 00946 Matrix<T> &mat; 00947 public: 00950 Matrix<T> operator*(const Matrix<T> &o) const { 00951 mat.checkDimensionMatch(o); 00952 Matrix<T> result(o.data->dim); 00953 MatrixCWiseMult::evaluate(mat, o, result); 00954 return result; 00955 } 00958 typename Matrix<T>::CWise& operator*=(const Matrix<T> &o){ 00959 mat.checkDimensionMatch(o); 00960 mat.data.detach(); 00961 MatrixCWiseMult::evaluate(mat, o, mat); 00962 return *this; 00963 } 00966 Matrix<T> operator/(const Matrix<T> &o) const { 00967 mat.checkDimensionMatch(o); 00968 Matrix<T> result(o.data->dim); 00969 MatrixCWiseDiv::evaluate(mat, o, result); 00970 return result; 00971 } 00974 typename Matrix<T>::CWise& operator/=(const Matrix<T> &o){ 00975 mat.checkDimensionMatch(o); 00976 mat.data.detach(); 00977 MatrixCWiseDiv::evaluate(mat, o, mat); 00978 return *this; 00979 } 00982 Matrix<T> max(const Matrix<T> &o) const { 00983 mat.checkDimensionMatch(o); 00984 Matrix<T> result(mat.data->dim); 00985 MatrixCWiseMax::evaluate(mat, o, result); 00986 return result; 00987 } 00990 Matrix<T> min(const Matrix<T> &o) const { 00991 mat.checkDimensionMatch(o); 00992 Matrix<T> result(mat.data->dim); 00993 MatrixCWiseMin::evaluate(mat, o, result); 00994 return result; 00995 } 00996 operator const Matrix<T>(){return mat;} 00997 }; 00998 01002 CWise cwise(){ 01003 return CWise(*this); 01004 } 01005 const CWise cwise() const { 01006 return CWise(*this); 01007 } 01008 01009 // ============================================================================= 01010 // Matrix unary operations 01011 // ============================================================================= 01012 private: 01013 #ifndef NEMO_PARSED_BY_DOXYGEN 01014 class MatrixMinus{ 01015 public: 01016 static void evaluate(const Matrix<T> &v, Matrix<T> &result){ 01017 result.data->mat = - v.data->mat; 01018 } 01019 }; 01020 #endif 01021 01022 public: 01023 Matrix<T> operator-() const { 01024 Matrix<T> result(dimension()); 01025 MatrixMinus::evaluate(*this, result); 01026 return result; 01027 } 01028 Matrix<T> operator+() const{ 01029 return *this; 01030 } 01031 01032 // ============================================================================= 01033 // Matrix-Vector operations 01034 // ============================================================================= 01035 private: 01036 #ifndef NEMO_PARSED_BY_DOXYGEN 01037 class MatrixMultVector{ 01038 public: 01039 static void evaluate(const Matrix<T> &m, const MathVector<T> &v, MathVector<T> &result){ 01040 // FIXME performance: no index checking!!! 01041 for(unsigned int row = 0; row < m.rows(); row++){ 01042 result[row] = 0.0; 01043 for(unsigned int col = 0; col < m.cols(); col++){ 01044 result[row] += m(row,col) * v[col]; 01045 } 01046 } 01047 } 01048 }; 01049 #endif 01050 01051 public: 01055 MathVector<T> operator*(const MathVector<T> &r) const{ 01056 checkColumnDimensionMatch(r); 01057 MathVector<T> result(dim(rows())); 01058 MatrixMultVector::evaluate(*this, r, result); 01059 return result; 01060 } 01061 01062 // ============================================================================= 01063 // Matrix-Scalar operations 01064 // ============================================================================= 01065 private: 01066 #ifndef NEMO_PARSED_BY_DOXYGEN 01067 class MatrixMultScalar{ 01068 public: 01069 static void evaluate(const Matrix<T> &m, const T &s, Matrix<T> &result){ 01070 for(unsigned int row = 0; row < m.rows(); row++){ 01071 for(unsigned int col = 0; col < m.cols(); col++){ 01072 // FIXME performance: no index checking!!! 01073 result(row,col) = m(row,col) * s; 01074 } 01075 } 01076 } 01077 }; 01078 #endif 01079 01080 public: 01084 Matrix<T> operator*(const T &r) const{ 01085 Matrix<T> result(this->dimension()); 01086 MatrixMultScalar::evaluate(*this, r, result); 01087 return result; 01088 } 01092 friend Matrix<T> operator*(const T &s, const Matrix<T> &m){ 01093 Matrix<T> result(m.dimension()); 01094 MatrixMultScalar::evaluate(m, s, result); 01095 return result; 01096 } 01100 Matrix<T>& operator*=(const T &r){ 01101 data.detach(); 01102 MatrixMultScalar::evaluate(*this, r, *this); 01103 return *this; 01104 } 01105 01106 private: 01107 01108 inline void checkColumnDimensionMatch(const Matrix<T> &r) const{ 01109 if(cols() != r.cols()){ 01110 _NEMO_ARGUMENT_ERROR("Matrix: dimension mismatch"); 01111 } 01112 } 01113 inline void checkColumnDimensionMatch(const MathVector<T> &v) const{ 01114 if(cols() != v.dimension()){ 01115 _NEMO_ARGUMENT_ERROR("Matrix/Vector: dimension mismatch"); 01116 } 01117 } 01118 inline void checkRowDimensionMatch(const Matrix<T> &r) const{ 01119 if(rows() != r.rows()){ 01120 _NEMO_ARGUMENT_ERROR("Matrix: dimension mismatch"); 01121 } 01122 } 01123 01124 /*inline void checkIndex(const unsigned int rowIndex, const unsigned int colIndex) const{ 01125 if(rows() <= rowIndex){ 01126 _NEMO_ARGUMENT_ERROR("Matrix: row-index out of range"); 01127 } 01128 if(cols() <= colIndex){ 01129 _NEMO_ARGUMENT_ERROR("Matrix: column-index out of range"); 01130 } 01131 }*/ 01132 template <class D> 01133 inline void checkDimensionMatch(const Matrix<D> &r) const{ 01134 if(dimension() != r.dimension()){ 01135 std::stringstream msg; 01136 msg << "Matrix: dimension mismatch: argument should have " 01137 << "dimension (" << dimension().rows() << "x" << dimension().cols() << ") " 01138 << "but has (" << r.dimension().rows() << "x" << r.dimension().cols() << ")"; 01139 _NEMO_ARGUMENT_ERROR(msg.str()); 01140 } 01141 } 01142 inline void checkSubsumes(const unsigned int fromRow, const unsigned int fromCol, MatrixDimension d) const { 01143 // check index (allowing index 0 for 0-dimension) 01144 if( this->rows() < fromRow+d.rows() ){ 01145 _NEMO_ARGUMENT_ERROR(""); 01146 } 01147 if( this->cols() < fromCol+d.cols() ){ 01148 _NEMO_ARGUMENT_ERROR(""); 01149 } 01150 } 01151 01152 IntrusiveCowSupport<MatrixData > data; 01153 01154 }; 01155 01156 // ============================================================================= 01157 // Inversion rules 01158 // ============================================================================= 01159 01160 #define _nemo_matrix_inversion_(TYPE) \ 01161 template <> inline Matrix<TYPE> leftAddInverse<Matrix<TYPE> >(const Matrix<TYPE> &value){return -value;} \ 01162 template <> inline Matrix<TYPE> rightAddInverse<Matrix<TYPE> >(const Matrix<TYPE> &value){return -value;} \ 01163 template <> inline Matrix<TYPE> leftSubtractInverse<Matrix<TYPE> >(const Matrix<TYPE> &value){return value;} \ 01164 template <> inline Matrix<TYPE> rightSubtractInverse<Matrix<TYPE> >(const Matrix<TYPE> &value){return -value;} 01165 01166 _nemo_matrix_inversion_(float) 01167 _nemo_matrix_inversion_(double) 01168 _nemo_matrix_inversion_(long double) 01169 //_nemo_matrix_inversion_(short) 01170 _nemo_matrix_inversion_(int) 01171 //_nemo_matrix_inversion_(long int) 01172 01173 01174 01175 01176 template <class T> 01177 std::ostream& operator<<(std::ostream& out, const Matrix<T> &mat){ 01178 01179 MatrixDimension dim = mat.dimension(); 01180 01181 if(dim.rows() == 0 || dim.cols() == 0){ 01182 out << std::string("[]"); 01183 } 01184 /*else if(dim.rows() == 0){ 01185 out << std::string("[]"); 01186 } 01187 else if(dim.cols() == 0){ 01188 out << std::string("["); 01189 for(unsigned int r = 0; r < dim.rows()-1; r++){ 01190 out << std::string(";") << std::endl; 01191 } 01192 out << std::string(";]"); 01193 }*/ 01194 else{ 01195 out << std::string("["); 01196 for(unsigned int r=0; r<dim.rows(); r++){ 01197 if(r!=0){ 01198 out << std::string(" "); 01199 } 01200 01201 for(unsigned int c=0; c<dim.cols()-1; c++){ 01202 out << mat(r,c) << ",\t"; 01203 } 01204 out << mat(r, dim.cols()-1); 01205 01206 if(r < dim.rows()-1){ 01207 out << ";" << std::endl; 01208 } 01209 } 01210 out << std::string("]"); 01211 } 01212 01213 return out; 01214 } 01215 01216 01217 typedef Matrix<double> RealMatrix; 01218 01219 01220 } // end namespace 01221 01222 01223 #endif 01224
Generated on Mon Feb 25 12:49:59 2013 for NemoMath by
![doxygen](doxygen.png)