Changeset 420 for trunk/lib/gslapi/matrix.h
 Timestamp:
 Dec 2, 2005, 1:15:50 AM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/lib/gslapi/matrix.h
r390 r420 1 1 // $Id$ 2 2 3 #ifndef _theplu_gslapi_matrix_ 4 #define _theplu_gslapi_matrix_ 3 #ifndef _theplu_gslapi_matrix_ 4 #define _theplu_gslapi_matrix_ 5 5 6 6 #include <c++_tools/gslapi/vector.h> 7 #include <c++_tools/utility/Exception.h> 7 8 8 9 #include <gsl/gsl_matrix.h> 9 10 #include <iostream> 10 11 11 12 12 namespace theplu { … … 20 20 /// well in the future. 21 21 /// 22 class matrix 23 { 24 public: 22 /// \par[File streams] Reading and writing vectors to file streams 23 /// are of course supported. These are implemented without using GSL 24 /// functionality, and thus binary read and write to streams are not 25 /// supported. 26 /// 27 /// \par[Matrix views] GSL matrix views are supported and these are 28 /// disguised as ordinary gslapi::matrix'es. A support function is 29 /// added, gslapi::matrix::isview(), that can be used to check if a 30 /// matrix object is a view. Note that view matricess do not own the 31 /// undelying data, and a view is not valid if the matrix owning the 32 /// data is deallocated. 33 /// 34 /// @note All GSL matrix related functions are not implement but 35 /// most functionality defined for GSL matrices can be achieved with 36 /// this interface class. Most notable GSL functionality not 37 /// supported are no binary file support and views on arrays, 38 /// gslapi::vectors, gsl_vectors, diagonals, subdiagonals, and 39 /// superdiagonals. If there is an interest from the user community, 40 /// the omitted functionality could be included. 41 /// 42 /// @todo Maybe it would be smart to create temporary objects need 43 /// for BLAS in constructors. As it is now temporary objects are 44 /// called before BLAS functionality iss used, cf. const matrix& 45 /// matrix::operator*=(const matrix& other) 46 /// 47 class matrix 48 { 49 public: 25 50 26 51 /// 27 52 /// The default constructor. 28 /// 29 matrix(void); 53 /// 54 /// This contructor does not initialize underlying (essential) 55 /// structures. 56 /// 57 inline matrix(void) : m_(NULL), view_(NULL) {} 30 58 31 59 /// 32 60 /// Constructor. Allocates memory space for \a r times \a c 33 61 /// elements, and sets all elements to \a init_value. 34 /// 35 matrix(const size_t& r, const size_t& c, const double init_value=0); 36 37 /// 62 /// 63 inline matrix(const size_t& r, const size_t& c, double init_value=0) 64 : view_(NULL) { m_ = gsl_matrix_alloc(r,c); set_all(init_value); } 65 66 /// 38 67 /// The copy constructor. 39 /// 40 matrix(const matrix&); 68 /// 69 /// @note If the object to be copied is a matrix view, the values 70 /// of the view will be copied, i.e. the view is not copied. 71 /// 72 inline matrix(const matrix& o) : view_(NULL) {m_=o.create_gsl_matrix_copy();} 73 74 /// 75 /// The matrix view constructor. 76 /// 77 /// Create a view of matrix \a m, with starting row \a offset_row, 78 /// starting column \a offset_column, row size \a n_row, and 79 /// column size \a n_column. 80 /// 81 /// A matrix view can be used as any matrix with the difference 82 /// that changes made to the view will also change the object that 83 /// is viewed. Also, using the copy constructor will create a new 84 /// matrix object that is a copy of whatever is viewed. If a copy 85 /// of the view is needed then you should use this constructor to 86 /// obtain a copy. 87 /// 88 /// @note If the object viewed by the view goes out of scope or is 89 /// deleted, the view becomes invalid and the result of further 90 /// use is undefined. 91 /// 92 matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row, 93 size_t n_column); 41 94 42 95 /// … … 46 99 /// column elements in a row must be separated with white space 47 100 /// characters except the new line (\\n) character. Rows, and 48 /// columns, are read conse quently and only rectangular matrices101 /// columns, are read consecutively and only rectangular matrices 49 102 /// are supported. Empty lines are ignored. End of input is at end 50 103 /// of file marker. 51 104 /// 52 explicit matrix(std::istream &) ;53 54 105 explicit matrix(std::istream &) throw (utility::IO_error,std::exception); 106 107 /// 55 108 /// The destructor. 56 109 /// 57 ~matrix(void); 58 59 /// 110 ~matrix(void); 111 112 /// 113 /// Elementwise addition of the elements of matrix \a b to the 114 /// elements of the calling matrix ,\f$a_{ij} = a_{ij} + b_{ij} \; 115 /// \forall i,j\f$. The result is stored into the calling matrix. 116 /// 117 /// @return Whatever GSL returns. 118 /// 119 inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); } 120 121 /// 122 /// Add the scalar value \a d to the elements of the calling 123 /// matrix, \f$a_{ij} = a_{ij} + d \; \forall i,j\f$. The result 124 /// is stored into the calling matrix. 125 /// 126 /// @return Whatever GSL returns. 127 /// 128 inline int 129 add_constant(const double d) { return gsl_matrix_add_constant(m_,d); } 130 131 /// 60 132 /// @return The number of columns in the matrix. 61 133 /// 62 inline size_t columns(void) const { return m_>size2; } 63 64 /// 65 /// Swap column \a i with \a j. 66 /// 67 inline void column_swap(const size_t& i,const size_t& j) 68 { gsl_matrix_swap_columns(m_,i,j); } 69 70 /// 71 /// @return true if each element deviates less or equal than \a 72 /// d. If any matrix contain a nan false is always returned. 73 /// 74 bool equal(const matrix&, const double d=0) const; 134 inline size_t columns(void) const { return m_>size2; } 135 136 /// 137 /// Elementwise division of the elemnts of the calling matrix by 138 /// the elements of matrix \a b, \f$a_{ij} = a_{ij} / b_{ij} \; 139 /// \forall i,j\f$. The result is stored into the calling matrix. 140 /// 141 /// @return Whatever GSL returns. 142 /// 143 inline int 144 div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); } 145 146 /// 147 /// Check whether matrices are equal within a user defined 148 /// precision, set by \a precision. 149 /// 150 /// @return True if each element deviates less or equal than \a 151 /// d. If any matrix contain a NaN, false is always returned. 152 /// 153 bool equal(const matrix&, const double precision=0) const; 154 155 /// 156 /// @return A const pointer to the internal GSL matrix. 157 /// 158 inline const gsl_matrix* gsl_matrix_p(void) const { return m_; } 159 160 /// 161 /// @return A pointer to the internal GSL matrix. 162 /// 163 // Jari, is this needed? 164 inline gsl_matrix* gsl_matrix_p(void) { return m_; }; 165 166 /// 167 /// @return True if all elements in the matrix is zero, false 168 /// othwerwise; 169 /// 170 inline bool isnull(void) const { return gsl_matrix_isnull(m_); } 171 172 /// 173 /// Check if the matrix object is a view (submatrix) to another 174 /// matrix. 175 /// 176 /// @return True if the object is a view, false othwerwise. 177 /// 178 inline bool isview(void) const { return view_; } 179 180 /// 181 /// @return The maximum value of the matrix. 182 /// 183 inline double max(void) const { return gsl_matrix_max(m_); } 184 185 /// 186 /// @return The index to the element with the maximum value of the 187 /// matrix. The lowest index has precedence (searching in 188 /// rowmajor order). 189 /// 190 inline std::pair<size_t,size_t> max_index(void) const; 191 192 /// 193 /// @return The minimum value of the matrix. 194 /// 195 inline double min(void) const { return gsl_matrix_min(m_); } 196 197 /// 198 /// @return The index to the element with the minimum value of the 199 /// matrix. The lowest index has precedence (searching in 200 /// rowmajor order). 201 /// 202 inline std::pair<size_t,size_t> min_index(void) const; 203 204 /// 205 /// @return The indecies to the element with the minimum and 206 /// maximum values of the matrix, respectively. The lowest index 207 /// has precedence (searching in rowmajor order). 208 /// 209 inline void minmax_index(std::pair<size_t,size_t>& min, 210 std::pair<size_t,size_t>& max) const 211 { gsl_matrix_minmax_index(m_,&min.first,&min.second, 212 &max.first,&max.second); } 213 214 /// 215 /// @return The minimum and maximum values of the matrix, as the 216 /// \a first and \a second member of the returned \a pair, 217 /// respectively. 218 /// 219 std::pair<double,double> minmax(void) const; 220 221 /// 222 /// Multiply the elements of matrix \a b with the elements of the 223 /// calling matrix ,\f$a_{ij} = a_{ij} * b_{ij} \; \forall 224 /// i,j\f$. The result is stored into the calling matrix. 225 /// 226 /// @return Whatever GSL returns. 227 /// 228 inline int 229 mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); } 230 231 /// 232 /// @return The number of rows in the matrix. 233 /// 234 inline size_t rows(void) const { return m_>size1; } 235 236 /// 237 /// Multiply the elements of the calling matrix with a scalar \a 238 /// d, \f$a_{ij} = d * a_{ij} \; \forall i,j\f$. The result is 239 /// stored into the calling matrix. 240 /// 241 /// @return Whatever GSL returns. 242 /// 243 inline int scale(const double d) { return gsl_matrix_scale(m_,d); } 244 245 /// 246 /// Set element values to values in \a mat. This function is 247 /// needed for assignment of viewed elements. 248 /// 249 /// @return Whatever GSL returns. 250 /// 251 /// @note No check on size is done. 252 /// 253 /// @see const matrix& operator=(const matrix&) 254 /// 255 inline int set(const matrix& mat) { return gsl_matrix_memcpy(m_,mat.m_); } 256 257 /// 258 /// Set all elements to \a value. 259 /// 260 inline void set_all(const double value) { gsl_matrix_set_all(m_,value); } 261 262 /// 263 /// Set \a column values to values in \a vec. 264 /// 265 /// @return Whatever GSL returns. 266 /// 267 /// @note No check on size is done. 268 /// 269 inline int set_column(const size_t column, const vector& vec) 270 { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); } 271 272 /// 273 /// Set \a row values to values in \a vec. 274 /// 275 /// @return Whatever GSL returns. 276 /// 277 /// @note No check on size is done. 278 /// 279 inline int set_row(const size_t row, const vector& vec) 280 { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); } 281 282 /// 283 /// Subtract the elements of matrix \a b from the elements of the 284 /// calling matrix ,\f$a_{ij} = a_{ij}  b_{ij} \; \forall 285 /// i,j\f$. The result is stored into the calling matrix. 286 /// 287 /// @return Whatever GSL returns. 288 /// 289 inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); } 290 291 /// 292 /// Exchange the elements of the this and \a other by copying. The 293 /// two matrices must have the same size. 294 /// 295 /// @return Whatever GSL returns. 296 /// 297 inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); } 298 299 /// 300 /// Swap columns \a i and \a j. 301 /// 302 inline int swap_columns(const size_t i,const size_t j) 303 { return gsl_matrix_swap_columns(m_,i,j); } 304 305 /// 306 /// Swap row \a i and column \a j. 307 /// 308 inline int swap_rowcol(const size_t i,const size_t j) 309 { return gsl_matrix_swap_rowcol(m_,i,j); } 310 311 /// 312 /// Swap rows \a i and \a j. 313 /// 314 inline int swap_rows(const size_t i, const size_t j) 315 { return gsl_matrix_swap_rows(m_,i,j); } 316 317 /// 318 /// Transpose the matrix. 319 /// 320 void transpose(void); 321 322 /// 323 /// @return Reference to the element position (\a row, \a column). 324 /// 325 inline double& operator()(size_t row,size_t column) 326 { return (*gsl_matrix_ptr(m_,row,column)); } 327 328 /// 329 /// @return Const reference to the element position (\a row, \a 330 /// column). 331 /// 332 inline const double& operator()(size_t row,size_t column) const 333 { return (*gsl_matrix_const_ptr(m_,row,column)); } 334 335 /// 336 /// Matrixvector multiplication. 337 /// 338 /// @return The resulting vector. 339 /// 340 // Jari, where should this go? 341 const vector operator*(const vector&) const; 342 343 /// 344 /// Comparison operator. 345 /// 346 /// @return True if all elements are equal otherwise False. 347 /// 348 /// @see equal 349 /// 350 inline bool operator==(const matrix& other) const { return equal(other); } 351 352 /// 353 /// Comparison operator. 354 /// 355 /// @return False if all elements are equal otherwise True. 356 /// 357 /// @see equal 358 /// 359 inline bool operator!=(const matrix& other) const { return !equal(other); } 360 361 /// 362 /// The assignment operator. There is no requirements on 363 /// dimensions, i.e. the matrix is remapped in memory if 364 /// necessary. This implies that in general views cannot be 365 /// assigned using this operator. Views will be mutated into 366 /// normal matrices. The only exception to this behaviour on views 367 /// is when selfassignemnt is done, since selfassignment is 368 /// ignored. 369 /// 370 /// @return A const reference to the resulting matrix. 371 /// 372 /// @see int set(const matrix&) 373 /// 374 const matrix& operator=(const matrix& other); 375 376 /// 377 /// Add and assign operator. 378 /// 379 inline const matrix& operator+=(const matrix& m) { add(m); return *this; } 380 381 /// 382 /// Add and assign operator. 383 /// 384 inline const matrix& 385 operator+=(const double d) { add_constant(d); return *this; } 386 387 /// 388 /// Subtract and assign operator. 389 /// 390 inline const matrix& operator=(const matrix& m) { sub(m); return *this; } 391 392 /// 393 /// Multiply and assigment operator. 394 /// 395 /// @return Const reference to the resulting matrix. 396 /// 397 const matrix& operator*=(const matrix&); 398 399 /// 400 /// Multiply and assign operator. 401 /// 402 inline const matrix& operator*=(const double d) { scale(d); return *this; } 403 404 405 private: 75 406 76 407 /// … … 82 413 /// @return A pointer to a copy of the internal GSL matrix. 83 414 /// 84 gsl_matrix* gsl_matrix_copy(void) const; 85 86 /// 87 /// @return A const pointer to the internal GSL matrix. 88 /// 89 inline const gsl_matrix* gsl_matrix_pointer(void) const { return m_; } 90 91 /// 92 /// @return A pointer to the internal GSL matrix. 93 /// 94 inline gsl_matrix* gsl_matrix_pointer(void) { return m_; }; 95 96 /// 97 /// This function multiplies all elements between two matrices and 98 /// the result is stored back into the calling object, \f$a_{ij} = 99 /// a_{ij} * b_{ij} \; \forall i,j\f$. 100 /// 101 inline void mul_elements(const matrix& b) 102 { gsl_matrix_mul_elements(m_,b.m_); } 103 104 /// 105 /// Swap row \a i with \a j. 106 /// 107 inline void row_swap( const size_t& i, const size_t& j) 108 { gsl_matrix_swap_rows(m_,i,j); } 109 110 /// 111 /// @return The number of rows in the matrix. 112 /// 113 inline size_t rows(void) const { return m_>size1; } 114 115 /// 116 /// Set all elements to \a value. 117 /// 118 inline void set_all(const double value) { gsl_matrix_set_all(m_,value); } 119 120 /// 121 /// Set column \a i to \a vec. 122 /// 123 inline void set_column(const size_t i, const vector vec) 124 {gsl_matrix_set_col(m_, i, vec.gsl_vector_pointer());} 125 126 /// 127 /// Set row \a i to \a vec. 128 /// 129 inline void set_row(const size_t i, const vector vec) 130 {gsl_matrix_set_row(m_, i, vec.gsl_vector_pointer());} 131 132 /// 133 /// Set column \a i to \a vec 134 /// 135 //void matrix set_column(const size_t i, const vector vec) const; 136 137 /// 138 /// Transpose the matrix. 139 /// 140 /// @return Transpose of the matrix. 141 /// 142 matrix transpose(void) const; 143 144 /// 145 /// @return Reference to the element position (\a row, \a column). 146 /// 147 inline double& operator()(size_t row,size_t column) 148 { return (*gsl_matrix_ptr(m_,row,column)); } 149 150 /// 151 /// @return Const reference to the element position (\a row, \a 152 /// column). 153 /// 154 inline const double& operator()(size_t row,size_t column) const 155 { return (*gsl_matrix_const_ptr(m_,row,column)); } 156 157 /// 158 /// @todo to be properly implemented. Should return a reference to a vector 159 /// object. Underlying GSL supports only GSL_vector views, thus 160 /// some redesign of the vector class is needed. 161 /// 162 /// inline vector& operator[](size_t i); 163 /// So for now, stupid copying is done, and actually a matrix 164 /// row is returned using this function. 165 vector operator[](size_t row) const; 166 167 /// 168 /// @todo to be properly implemented. Should return a reference to 169 /// a vector object (matrix row). Underlying GSL supports only GSL_vector 170 /// views, thus some redesign of the vector class is needed. 171 /// 172 /// inline const vector& operator[](size_t i) const; 173 /// 174 /// So for now, stupid copying is done, and actually a matrix 175 /// column is returned using this function. 176 vector TEMP_col_return(size_t column) const; 177 178 /// 179 /// Matrix multiplication. 180 /// 181 /// @return The resulting matrix. 182 /// 183 matrix operator*(const matrix&) const; 184 185 /// 186 /// Matrixvector multiplication. 187 /// 188 /// @return The resulting vector. 189 /// 190 vector operator*(const vector&) const; 191 192 /// 193 /// Matrix addition. 194 /// 195 /// @return The resulting matrix. 196 /// 197 matrix operator+(const matrix&) const; 198 199 /// 200 /// @todo Matrix addition. 201 /// 202 /// @return The resulting matrix. 203 /// 204 matrix operator+=(const matrix&); 205 206 /// 207 /// Matrix subtraction. 208 /// 209 /// @return The resulting matrix. 210 /// 211 matrix operator(const matrix&) const; 212 213 /// 214 /// @todo Matrix subtraction. 215 /// 216 /// @return The resulting matrix. 217 /// 218 matrix operator=(const matrix&) const; 219 220 /// 221 /// Comparison operator. 222 /// 223 /// @return True if all elements are equal otherwise False. 224 /// 225 /// @see equal 226 /// 227 inline bool operator==( const matrix& other ) const { return equal(other);} 228 229 /// 230 /// Comparison operator. 231 /// 232 /// @return False if all elements are equal otherwise True. 233 /// 234 /// @see equal 235 /// 236 inline bool operator!=( const matrix& other ) const {return !equal(other);} 237 238 /// 239 /// The assignment operator. There is no requirements on 240 /// dimensions. 241 /// 242 matrix& operator=(const matrix& other); 243 244 /// 245 /// Subtract and assign operator. 246 /// 247 inline matrix& operator=(const matrix& m) 248 { gsl_matrix_sub(m_,m.m_); return *this; } 249 250 /// 251 /// Multiply and assign operator. 252 /// 253 inline matrix& operator*=(const double d) 254 { gsl_matrix_scale(m_,d); return *this; } 255 256 257 private: 258 259 gsl_matrix* m_; 260 415 gsl_matrix* create_gsl_matrix_copy(void) const; 416 417 gsl_matrix* m_; 418 gsl_matrix_view* view_; 261 419 }; 262 420 263 421 /// 264 422 /// The output operator for the matrix class. 265 /// 266 std::ostream& operator<< (std::ostream& s, const matrix&); 267 423 /// 424 std::ostream& operator<< (std::ostream& s, const matrix&); 268 425 269 426 270 427 }} // of namespace gslapi and namespace theplu 271 428 272 #endif 429 #endif
Note: See TracChangeset
for help on using the changeset viewer.