00001 /*************************************************************************** 00002 ** 00003 ** This file is part of QGpCoreTools. 00004 ** 00005 ** This library is free software; you can redistribute it and/or 00006 ** modify it under the terms of the GNU Lesser General Public 00007 ** License as published by the Free Software Foundation; either 00008 ** version 2.1 of the License, or (at your option) any later version. 00009 ** 00010 ** This file is distributed in the hope that it will be useful, but WITHOUT 00011 ** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 ** FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 00013 ** License for more details. 00014 ** 00015 ** You should have received a copy of the GNU Lesser General Public 00016 ** License along with this library; if not, write to the Free Software 00017 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00018 ** 00019 ** See http://www.geopsy.org for more information. 00020 ** 00021 ** Created : 2004-05-17 00022 ** Authors : 00023 ** Marc Wathelet 00024 ** Marc Wathelet (ULg, Liège, Belgium) 00025 ** Marc Wathelet (LGIT, Grenoble, France) 00026 ** 00027 ***************************************************************************/ 00028 00029 #ifndef MATRIX_H 00030 #define MATRIX_H 00031 00032 00033 #include "MatrixData.h" 00034 #include "Point.h" 00035 #include "Complex.h" 00036 00037 namespace QGpCoreTools { 00038 00039 template <typename T> class MatrixIterator; 00040 template <typename T> class MutableMatrixIterator; 00041 template <typename T> class MatrixRowIterator; 00042 template <typename T> class MutableMatrixRowIterator; 00043 template <typename T> class MatrixColumnIterator; 00044 template <typename T> class MutableMatrixColumnIterator; 00045 00046 template <typename T> class QGPCORETOOLS_EXPORT Matrix 00047 { 00048 public: 00049 Matrix() {_d=new MatrixData<T>;} 00050 Matrix(int ndim) {_d=new MatrixData<T>; _d->resize(ndim);} 00051 Matrix(int nrow, int ncol) {_d=new MatrixData<T>; _d->resize(nrow, ncol);} 00052 Matrix(const Matrix<T>& o) : _d(o._d) {} 00053 00054 Matrix<T> operator=(const Matrix<T>& m) {_d=m._d; return *this;} 00055 bool operator==(const Matrix<T>& m) const {return _d->operator==(*m._d);} 00056 bool operator!=(const Matrix<T>& m) const {return !_d->operator==(*m._d);} 00057 Matrix<T> operator+(const Matrix<T>& m) const {return _d->operator+(*m._d);} 00058 void operator+=(const Matrix<T>& m); 00059 Matrix<T> operator-(const Matrix<T>& m) const {return _d->operator-(*m._d);} 00060 void operator-=(const Matrix<T>& m); 00061 Matrix<T> operator*(const Matrix<T>& m) const {return _d->operator*(*m._d);} 00062 void operator*=(const Matrix<T>& m); 00063 void zero() {_d->zero();} 00064 void identity() {_d->identity();} 00065 void fill(const T& v) {_d->fill(v);} 00066 T& at(int row, int col) {return _d->at(row,col);} 00067 const T& at(int row, int col) const {return _d->at(row,col);} 00068 const T * data() const {return _d->data();} 00069 T * data() {return _d->data();} 00070 QVector<T> rowAt(int row) const; 00071 QVector<T> columnAt(int col) const; 00072 Matrix<T> transposed() const {return _d->transposed();} 00073 void invert() {_d->invert();} 00074 void transpose() { *this=transposed();} 00075 void mergeRow(const Matrix<T>& row1, const Matrix<T>& row2) {_d->mergeRow(*row1._d, *row2._d);} 00076 void mergeColumn(const Matrix<T>& col1, const Matrix<T>& col2) {_d->mergeColumn(*col1._d, *col2._d);} 00077 void copyAt(int row, int col, const Matrix<T>& m) {_d->copyAt(row, col, *m._d);} 00078 Matrix<T> subMatrix(int rowStart, int colStart, int rowEnd, int colEnd) const {return _d->subMatrix(rowStart, colStart, rowEnd, colEnd);} 00079 Matrix<T> sortedRows(const PermutationVector& p) {return _d->sortedRows(p);} 00080 Matrix<T> sortedColumns(const PermutationVector& p) {return _d->sortedColumns(p);} 00081 00082 int columnCount() const {return _d->columnCount();} 00083 int rowCount() const {return _d->rowCount();} 00084 int nextNullRow(int startRow) const {return _d->nextNullRow(startRow);} 00085 void resize(int ndim) {_d->resize(ndim);} 00086 void resize(int nrow, int ncol) {_d->resize(nrow, ncol);} 00087 00088 QString toUserString(int precision=6, char format='g') const {return _d->toUserString(precision, format);} 00089 QString toString() const {return _d->toString();} 00090 bool fromString(const StringSection& s) {return _d->fromString(s);} 00091 protected: 00092 void swapRowColumn() {_d->swapRowColumn();} 00093 00094 friend class MatrixData<T>; 00095 friend class MatrixIterator<T>; 00096 friend class MutableMatrixIterator<T>; 00097 friend class MatrixRowIterator<T>; 00098 friend class MutableMatrixRowIterator<T>; 00099 friend class MatrixColumnIterator<T>; 00100 friend class MutableMatrixColumnIterator<T>; 00101 00102 QSharedDataPointer< MatrixData<T> > _d; 00103 }; 00104 00105 template <typename T> Matrix<T> MatrixData<T>::operator+(const MatrixData<T>& m) const 00106 { 00107 ASSERT(_nrow==m._nrow && _ncol==m._ncol); 00108 Matrix<T> r(_nrow, _ncol); 00109 T * rValues=r._d->values(); 00110 for(int i=_nrow*_ncol-1; i>=0; i--) { 00111 rValues[i]=_values[i]+m._values[i]; 00112 } 00113 return r; 00114 } 00115 00116 template <typename T> Matrix<T> MatrixData<T>::operator-(const MatrixData<T>& m) const 00117 { 00118 ASSERT(_nrow==m._nrow && _ncol==m._ncol); 00119 Matrix<T> r(_nrow, _ncol); 00120 T * rValues=r._d->values(); 00121 for(int i=_nrow*_ncol-1; i>=0; i--) { 00122 rValues[i]=_values[i]-m._values[i]; 00123 } 00124 return r; 00125 } 00126 00127 template <typename T> Matrix<T> MatrixData<T>::operator*(const MatrixData<T>& m) const 00128 { 00129 ASSERT(_ncol==m._nrow); 00130 Matrix<T> r(_nrow, m._ncol); 00131 r.zero(); 00132 T * rValues=r._d->_values; 00133 const T * thisValues=_values; 00134 const T * mValues=m._values; 00135 if(MatrixMultiply::isValid(_nrow, _ncol, m._ncol)) { 00136 const MatrixMultiply& multiplier=*MatrixMultiply::begin(_nrow, _ncol, m._ncol); 00137 int n=multiplier.indexCount(); 00138 for(int i=0;i<n;i++) { 00139 const MatrixMultiply::IndexMap& index=multiplier.indexMap(i); 00140 rValues[index.resultIndex()]+=thisValues[index.factor1Index()]*mValues[index.factor2Index()]; 00141 } 00142 MatrixMultiply::end(&multiplier); 00143 } else { 00144 // Slower method because cache would be too big for matrix>256*256 00145 for(int i=0; i<_nrow; i++) { 00146 for(int j=0; j<m._ncol; j++) { 00147 int resultIndex=i+j*_nrow; 00148 int mIndex=j*_ncol; 00149 for(int k=0; k<_ncol; k++) { 00150 rValues[resultIndex]+=thisValues[i+k*_nrow]*mValues[k+mIndex]; 00151 } 00152 } 00153 } 00154 } 00155 return r; 00156 } 00157 00158 template <typename T> Matrix<T> MatrixData<T>::transposed() const 00159 { 00160 Matrix<T> t(_ncol, _nrow); 00161 T * pt=t._d->_values; 00162 for(int i=0; i<_nrow; i++) { 00163 for(int j=0; j<_ncol; j++) { 00164 (*pt++)=_values[j*_nrow+i]; 00165 } 00166 } 00167 return t; 00168 } 00169 00170 template <typename T> Matrix<T> MatrixData<T>::subMatrix(int rowStart, int colStart, int rowEnd, int colEnd) const 00171 { 00172 ASSERT(rowStart>=0 && rowStart<=rowEnd && rowEnd<_nrow); 00173 ASSERT(colStart>=0 && colStart<=colEnd && colEnd<_ncol); 00174 int mnrow=rowEnd-rowStart+1; 00175 int mncol=colEnd-colStart+1; 00176 Matrix<T> m(mnrow, mncol); 00177 int length=sizeof(T)*mnrow; 00178 T * values=m._d->_values; 00179 T * lastValue=values+mnrow*mncol; 00180 T * thisValues=_values+rowStart+colStart*_nrow; 00181 while(values<lastValue) { 00182 memcpy(values, thisValues, length); 00183 values+=mnrow; 00184 thisValues+=_nrow; 00185 } 00186 return m; 00187 } 00188 00189 template <typename T> Matrix<T> MatrixData<T>::sortedRows(const PermutationVector& p) 00190 { 00191 Matrix<T> m(_nrow, _ncol); 00192 for(int i=0; i<_ncol; i++) { 00193 int offset=i*_nrow; 00194 T * values=m._d->_values+offset; 00195 T * thisValues=_values+offset; 00196 for(int j=0; j<_nrow; j++) { 00197 values[p.newIndex(j)]=thisValues[j]; 00198 } 00199 } 00200 return m; 00201 } 00202 00203 template <typename T> Matrix<T> MatrixData<T>::sortedColumns(const PermutationVector& p) 00204 { 00205 Matrix<T> m(_nrow, _ncol); 00206 T * values=m._d->_values; 00207 for(int i=0; i<_ncol; i++) { 00208 memcpy(values+p.newIndex(i)*_nrow, _values+i*_nrow, sizeof(T)*_nrow); 00209 } 00210 return m; 00211 } 00212 00213 template <typename T> inline void Matrix<T>::operator*=(const Matrix<T>& m) 00214 { 00215 *this=_d->operator*( *m._d); 00216 } 00217 00218 template <typename T> inline void Matrix<T>::operator+=(const Matrix<T>& m) 00219 { 00220 *this=_d->operator+( *m._d); 00221 } 00222 00223 template <typename T> inline void Matrix<T>::operator-=(const Matrix<T>& m) 00224 { 00225 *this=_d->operator-( *m._d); 00226 } 00227 00228 template <typename T> QVector<T> Matrix<T>::rowAt(int row) const 00229 { 00230 QVector<T> r(columnCount()); 00231 MatrixColumnIterator<T> it(*this, row); 00232 int i=0; 00233 while(it.hasNext()) { 00234 r[i++]=*it.next(); 00235 } 00236 return r; 00237 } 00238 00239 template <typename T> QVector<T> Matrix<T>::columnAt(int col) const 00240 { 00241 QVector<T> r(rowCount()); 00242 MatrixColumnIterator<T> it(*this, col); 00243 int i=0; 00244 while(it.hasNext()) { 00245 r[i++]=*it.next(); 00246 } 00247 return r; 00248 } 00249 00250 00251 template <typename T> QDataStream& operator<<(QDataStream& s, const Matrix<T>& m) 00252 { 00253 s << m.rowCount() << m.columnCount(); 00254 for(int i=m.rowCount()*m.columnCount()-1; i>=0; i--) { 00255 s << m.data()[i]; 00256 } 00257 return s; 00258 } 00259 00260 template <typename T> QDataStream& operator>>(QDataStream& s, Matrix<T>& m) 00261 { 00262 int nrow, ncol; 00263 s >> nrow >> ncol; 00264 m.resize(nrow, ncol); 00265 for(int i=nrow*ncol-1; i>=0; i--) { 00266 s >> m.data()[i]; 00267 } 00268 return s; 00269 } 00270 00271 } // namespace QGpCoreTools 00272 00273 #endif // MATRIX_H
