Changeset 403
 Timestamp:
 Nov 27, 2005, 1:26:25 AM (17 years ago)
 Location:
 branches/better_matrix_class/lib/gslapi
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

branches/better_matrix_class/lib/gslapi/matrix.cc
r357 r403 53 53 std::vector<double> v; 54 54 for (nof_rows = 0; utility::read_to_double(is, v); nof_rows++) { 55 55 56 56 // Ignoring empty lines 57 57 if(!v.size()){ … … 59 59 continue; 60 60 } 61 61 62 62 if(nof_columns==0) 63 63 nof_columns=v.size(); 64 64 else if(v.size()!=nof_columns) { 65 65 std::cerr << "matrix::matrix(std::istream&) data file error: " 66 << "line" << nof_rows+1 << " has " << v.size() 66 << "line" << nof_rows+1 << " has " << v.size() 67 67 << " columns; expected " << nof_columns 68 68 << " columns" 69 << std::endl; 69 << std::endl; 70 70 exit(1); 71 71 } … … 77 77 // convert the data to a gsl matrix 78 78 m_ = gsl_matrix_alloc ( nof_rows, nof_columns ); 79 for(u_int i=0;i<nof_rows;i++) 80 for(u_int j=0;j<nof_columns;j++) 79 for(u_int i=0;i<nof_rows;i++) 80 for(u_int j=0;j<nof_columns;j++) 81 81 gsl_matrix_set( m_, i, j, data_matrix[i][j] ); 82 82 } … … 94 94 95 95 96 gsl_matrix* matrix::gsl_matrix_copy(void) const97 {98 gsl_matrix* m = gsl_matrix_alloc(rows(),columns());99 gsl_matrix_memcpy(m,m_);100 return m;101 }102 103 96 bool matrix::equal(const matrix& other, const double d) const 104 97 { 105 98 if (columns()!=other.columns()  rows()!=other.rows()) 106 99 return false; 107 for (size_t i=0; i<rows(); i++) 108 for (size_t j=0; j<columns(); j++) 100 for (size_t i=0; i<rows(); i++) 101 for (size_t j=0; j<columns(); j++) 109 102 // two last cond. are to check for nans 110 if (fabs( (*this)(i,j)other(i,j) ) > d  103 if (fabs( (*this)(i,j)other(i,j) ) > d  111 104 (*this)(i,j)!=(*this)(i,j)  other(i,j)!=other(i,j)) 112 105 return false; … … 114 107 return true; 115 108 } 109 110 111 112 gsl_matrix* matrix::gsl_matrix_copy(void) const 113 { 114 gsl_matrix* m = gsl_matrix_alloc(rows(),columns()); 115 gsl_matrix_memcpy(m,m_); // Jari, a GSL return value is ignored here 116 return m; 117 } 118 119 116 120 117 121 vector matrix::TEMP_col_return(size_t column) const … … 172 176 173 177 174 178 175 179 matrix& matrix::operator=( const matrix& other ) 176 180 { … … 179 183 gsl_matrix_free( m_ ); 180 184 m_ = other.gsl_matrix_copy(); 181 185 } 182 186 return *this; 183 187 } … … 185 189 186 190 191 // Jari, checkout GSL transpose support in GSL manual 8.4.9 187 192 matrix matrix::transpose(void) const 188 193 { … … 199 204 s.setf( ios::dec ); 200 205 s.precision(12); 201 for(size_t i=0, j=0; i<m.rows(); i++) 206 for(size_t i=0, j=0; i<m.rows(); i++) 202 207 for (j=0; j<m.columns(); j++) { 203 208 s << m(i,j); … … 207 212 s << "\n"; 208 213 } 209 214 return s; 210 215 } 211 216 
branches/better_matrix_class/lib/gslapi/matrix.h
r390 r403 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> … … 20 20 /// well in the future. 21 21 /// 22 class matrix 23 { 24 public: 22 /// @note All GSL matrix related functions are not implement but all 23 /// functionality, except binary file support, defined for GSL 24 /// matrices can be achieved with this interface class. 25 /// 26 class matrix 27 { 28 public: 25 29 26 30 /// 27 31 /// The default constructor. 28 29 32 /// 33 matrix(void); 30 34 31 35 /// 32 36 /// Constructor. Allocates memory space for \a r times \a c 33 37 /// elements, and sets all elements to \a init_value. 34 35 36 37 38 /// 39 matrix(const size_t& r, const size_t& c, const double init_value=0); 40 41 /// 38 42 /// The copy constructor. 39 /// 40 matrix(const matrix&); 43 /// 44 /// @note If the object to be copied is a matrix view, the values 45 /// of the view will be copied, i.e. the view is not copied. 46 /// 47 matrix(const matrix&); 41 48 42 49 /// … … 46 53 /// column elements in a row must be separated with white space 47 54 /// characters except the new line (\\n) character. Rows, and 48 /// columns, are read conse quently and only rectangular matrices55 /// columns, are read consecutively and only rectangular matrices 49 56 /// are supported. Empty lines are ignored. End of input is at end 50 57 /// of file marker. … … 52 59 explicit matrix(std::istream &); 53 60 54 61 /// 55 62 /// The destructor. 56 63 /// 57 ~matrix(void); 58 59 /// 64 ~matrix(void); 65 66 /// 67 /// Elementwise addition of the elements of matrix \a b to the 68 /// elements of the calling matrix ,\f$a_{ij} = a_{ij} + b_{ij} \; 69 /// \forall i,j\f$. The result is stored into the calling matrix. 70 /// 71 /// @return Whatever GSL returns. 72 /// 73 inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); } 74 75 /// 76 /// Add the scalar value \a d to the elements of the calling 77 /// matrix, \f$a_{ij} = a_{ij} + d \; \forall i,j\f$. The result 78 /// is stored into the calling matrix. 79 /// 80 /// @return Whatever GSL returns. 81 /// 82 inline int add_constant(const double d) 83 { return gsl_matrix_add_constant(m_,d); } 84 85 /// 60 86 /// @return The number of columns in the matrix. 61 87 /// 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); } 88 inline size_t columns(void) const { return m_>size2; } 89 90 /// 91 /// Elementwise division of the elemnts of the calling matrix by 92 /// the elements of matrix \a b, \f$a_{ij} = a_{ij} / b_{ij} \; 93 /// \forall i,j\f$. The result is stored into the calling matrix. 94 /// 95 /// @return Whatever GSL returns. 96 /// 97 inline int div_elements(const matrix& b) 98 { return gsl_matrix_div_elements(m_,b.m_); } 69 99 70 100 /// … … 75 105 76 106 /// 77 /// Create a new copy of the internal GSL matrix. 78 /// 79 /// Necessary memory for the new GSL matrix is allocated and the 80 /// caller is responsible for freeing the allocated memory. 81 /// 82 /// @return A pointer to a copy of the internal GSL matrix. 83 /// 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 /// 107 /// @return True if all elements in the matrix is zero, false 108 /// othwerwise; 109 /// 110 inline bool isnull(void) const { return gsl_matrix_isnull(m_); } 111 112 /// 113 /// @return The maximum value of the matrix. 114 /// 115 inline double max(void) const { return gsl_matrix_max(m_); } 116 117 /// 118 /// @return The index to the element with the maximum value of the 119 /// matrix. The lowest index has precedence (searching in 120 /// rowmajor order). 121 /// 122 inline std::pair<size_t,size_t> max_index(void) const; 123 124 /// 125 /// @return The minimum value of the matrix. 126 /// 127 inline double min(void) const { return gsl_matrix_min(m_); } 128 129 /// 130 /// @return The index to the element with the minimum value of the 131 /// matrix. The lowest index has precedence (searching in 132 /// rowmajor order). 133 /// 134 inline std::pair<size_t,size_t> min_index(void) const; 135 136 /// 137 /// @return The indecies to the element with the minimum and 138 /// maximum values of the matrix, respectively. The lowest index 139 /// has precedence (searching in rowmajor order). 140 /// 141 inline void minmax_index(std::pair<size_t,size_t>& min, 142 std::pair<size_t,size_t>& max) const 143 { gsl_matrix_minmax_index(m_,&min.first,&min.second, 144 &max.first,&max.second); } 145 146 /// 147 /// @return The minimum and maximum values of the matrix, as the 148 /// \a first and \a second member of the returned \a pair, 149 /// respectively. 150 /// 151 std::pair<double,double> minmax(void) const; 152 153 /// 154 /// Multiply the elements of matrix \a b with the elements of the 155 /// calling matrix ,\f$a_{ij} = a_{ij} * b_{ij} \; \forall 156 /// i,j\f$. The result is stored into the calling matrix. 157 /// 158 /// @return Whatever GSL returns. 159 /// 160 inline int mul_elements(const matrix& b) 161 { return gsl_matrix_mul_elements(m_,b.m_); } 162 163 /// 111 164 /// @return The number of rows in the matrix. 112 /// 113 inline size_t rows(void) const { return m_>size1; } 165 /// 166 inline size_t rows(void) const { return m_>size1; } 167 168 /// 169 /// Multiply the elements of the calling matrix with a scalar \a 170 /// d, \f$a_{ij} = d * a_{ij} \; \forall i,j\f$. The result is 171 /// stored into the calling matrix. 172 /// 173 /// @return Whatever GSL returns. 174 /// 175 inline int scale(const double d) { return gsl_matrix_scale(m_,d); } 114 176 115 177 /// 116 178 /// 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 /// 179 /// 180 inline void set_all(const double value) { gsl_matrix_set_all(m_,value); } 181 182 /// 183 /// Set \a column to \a vec. 184 /// 185 /// @return Whatever GSL returns. 186 /// 187 inline int set_column(const size_t column, const vector& vec) 188 { return gsl_matrix_set_col(m_, column, vec.gsl_vector_pointer()); } 189 190 /// 191 /// Set \a row to \a vec. 192 /// 193 /// @return Whatever GSL returns. 194 /// 195 inline int set_row(const size_t row, const vector& vec) 196 { return gsl_matrix_set_row(m_, row, vec.gsl_vector_pointer()); } 197 198 /// 199 /// Subtract the elements of matrix \a b from the elements of the 200 /// calling matrix ,\f$a_{ij} = a_{ij}  b_{ij} \; \forall 201 /// i,j\f$. The result is stored into the calling matrix. 202 /// 203 /// @return Whatever GSL returns. 204 /// 205 inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); } 206 207 /// 208 /// Exchange the elements of the this and \a other by copying. The 209 /// two matrices must have the same size. 210 /// 211 /// @return Whatever GSL returns. 212 /// 213 inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); } 214 215 /// 216 /// Swap columns \a i and \a j. 217 /// 218 inline int swap_columns(const size_t i,const size_t j) 219 { return gsl_matrix_swap_columns(m_,i,j); } 220 221 /// 222 /// Swap row \a i and column \a j. 223 /// 224 inline int swap_rowcol(const size_t i,const size_t j) 225 { return gsl_matrix_swap_rowcol(m_,i,j); } 226 227 /// 228 /// Swap rows \a i and \a j. 229 /// 230 inline int swap_rows(const size_t i, const size_t j) 231 { return gsl_matrix_swap_rows(m_,i,j); } 232 233 /// 138 234 /// Transpose the matrix. 139 235 /// 140 /// @return Transpose of the matrix. 141 /// 142 matrix transpose(void) const; 143 236 void transpose(void); 237 144 238 /// 145 239 /// @return Reference to the element position (\a row, \a column). 146 240 /// 147 241 inline double& operator()(size_t row,size_t column) 148 242 { return (*gsl_matrix_ptr(m_,row,column)); } … … 151 245 /// @return Const reference to the element position (\a row, \a 152 246 /// column). 153 247 /// 154 248 inline const double& operator()(size_t row,size_t column) const 155 249 { return (*gsl_matrix_const_ptr(m_,row,column)); } 156 250 157 251 /// 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; 252 /// Write me, should return a vector view. 253 /// 254 vector operator[](size_t row) const; 166 255 167 256 /// … … 169 258 /// a vector object (matrix row). Underlying GSL supports only GSL_vector 170 259 /// views, thus some redesign of the vector class is needed. 171 172 260 /// 261 /// inline const vector& operator[](size_t i) const; 173 262 /// 174 263 /// So for now, stupid copying is done, and actually a matrix 175 264 /// column is returned using this function. 265 /// 266 /// Jari, cf. the first two function in GSL manual section 8.4.8 176 267 vector TEMP_col_return(size_t column) const; 177 268 178 /// 179 /// Matrix multiplication. 180 /// 181 /// @return The resulting matrix. 182 /// 183 matrix operator*(const matrix&) const; 184 185 /// 269 /// 186 270 /// Matrixvector multiplication. 187 271 /// 188 272 /// @return The resulting vector. 189 273 /// 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; 274 // Jari, where should this go? 275 const vector operator*(const vector&) const; 219 276 220 277 /// … … 225 282 /// @see equal 226 283 /// 227 inline bool operator==( const matrix& other ) const { return equal(other);}284 inline bool operator==(const matrix& other) const { return equal(other); } 228 285 229 286 /// … … 234 291 /// @see equal 235 292 /// 236 inline bool operator!=( const matrix& other ) const {return !equal(other);}293 inline bool operator!=(const matrix& other) const { return !equal(other); } 237 294 238 295 /// … … 240 297 /// dimensions. 241 298 /// 242 matrix& operator=(const matrix& other); 299 const matrix& operator=(const matrix& other); 300 301 /// 302 /// Add and assign operator. 303 /// 304 inline const matrix& operator+=(const matrix& m) { add(m); return *this; } 305 306 /// 307 /// Add and assign operator. 308 /// 309 inline const matrix& operator+=(const double d) 310 { add_constant(d); return *this; } 243 311 244 312 /// 245 313 /// Subtract and assign operator. 246 314 /// 247 inline matrix& operator=(const matrix& m) 248 { gsl_matrix_sub(m_,m.m_); return *this; } 315 inline const matrix& operator=(const matrix& m) { sub(m); return *this; } 316 317 /// 318 /// Multiply and assigment operator. 319 /// 320 /// @return Const reference to the resulting matrix. 321 /// 322 const matrix& operator*=(const matrix&); 249 323 250 324 /// 251 325 /// Multiply and assign operator. 252 326 /// 253 inline matrix& operator*=(const double d) 254 { gsl_matrix_scale(m_,d); return *this; } 255 256 257 private: 258 259 gsl_matrix* m_; 327 inline const matrix& operator*=(const double d) { scale(d); return *this; } 328 329 330 private: 331 /// 332 /// Create a new copy of the internal GSL matrix. 333 /// 334 /// Necessary memory for the new GSL matrix is allocated and the 335 /// caller is responsible for freeing the allocated memory. 336 /// 337 /// @return A pointer to a copy of the internal GSL matrix. 338 /// 339 gsl_matrix* gsl_matrix_copy(void) const; 340 341 /// 342 /// @return A const pointer to the internal GSL matrix. 343 /// 344 // Jari, is this needed? 345 inline const gsl_matrix* gsl_matrix_pointer(void) const { return m_; } 346 347 /// 348 /// @return A pointer to the internal GSL matrix. 349 /// 350 // Jari, is this needed? 351 inline gsl_matrix* gsl_matrix_pointer(void) { return m_; }; 352 353 gsl_matrix* m_; 260 354 261 355 }; … … 263 357 /// 264 358 /// The output operator for the matrix class. 265 /// 266 std::ostream& operator<< (std::ostream& s, const matrix&); 267 268 359 /// 360 std::ostream& operator<< (std::ostream& s, const matrix&); 269 361 270 362 }} // of namespace gslapi and namespace theplu 271 363 272 #endif 364 #endif
Note: See TracChangeset
for help on using the changeset viewer.