Changeset 792 for trunk/yat/utility
 Timestamp:
 Mar 11, 2007, 1:14:24 AM (15 years ago)
 Location:
 trunk/yat/utility
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

trunk/yat/utility/PCA.cc
r789 r792 63 63 std::auto_ptr<SVD> pSVD( new SVD( A_center ) ); 64 64 pSVD>decompose(); 65 utility::matrix U = pSVD>U();66 utility::matrix V = pSVD>V();65 utility::matrix U(pSVD>U()); 66 utility::matrix V(pSVD>V()); 67 67 68 68 // Read the eigenvectors and eigenvalues 69 eigenvectors_ = U;69 eigenvectors_.clone(U); 70 70 eigenvectors_ .transpose(); 71 71 eigenvalues_.clone(pSVD>s()); … … 103 103 std::auto_ptr<SVD> pSVD( new SVD( A_center ) ); 104 104 pSVD>decompose(); 105 utility::matrix U = pSVD>U();106 utility::matrix V = pSVD>V();105 utility::matrix U(pSVD>U()); 106 utility::matrix V(pSVD>V()); 107 107 108 108 // Read the eigenvectors and eigenvalues 109 eigenvectors_ =V;109 eigenvectors_.clone(V); 110 110 eigenvectors_.transpose(); 111 111 eigenvalues_.clone(pSVD>s()); … … 197 197 { 198 198 size_t DIM = eigenvalues_.size(); 199 explained_intensity_ = utility::vector( DIM);199 explained_intensity_.clone(utility::vector( DIM )); 200 200 double sum = 0; 201 201 for( size_t i = 0; i < DIM; ++i ) 
trunk/yat/utility/matrix.cc
r788 r792 29 29 #include "utility.h" 30 30 31 #include <cassert> 31 32 #include <cmath> 32 33 #include <sstream> … … 41 42 42 43 matrix::matrix(void) 43 : blas_result_(NULL), m_(NULL), view_(NULL) 44 : blas_result_(NULL), m_(NULL), view_(NULL), view_const_(NULL), 45 proxy_m_(NULL) 44 46 { 45 47 } … … 47 49 48 50 matrix::matrix(const size_t& r, const size_t& c, double init_value) 49 : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL) 51 : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL), 52 view_const_(NULL), proxy_m_(m_) 50 53 { 51 54 if (!m_) … … 56 59 57 60 matrix::matrix(const matrix& o) 58 : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL) 61 : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL), 62 view_const_(NULL), proxy_m_(m_) 59 63 { 60 64 } … … 63 67 matrix::matrix(matrix& m, size_t offset_row, size_t offset_column, 64 68 size_t n_row, size_t n_column) 65 : blas_result_(NULL) 69 : blas_result_(NULL), view_const_(NULL) 66 70 { 67 71 view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_, … … 70 74 if (!view_) 71 75 throw utility::GSL_error("matrix::matrix failed to setup view"); 72 m_ = &(view_>matrix);76 proxy_m_ = m_ = &(view_>matrix); 73 77 } 74 78 … … 77 81 matrix::matrix(std::istream& is, char sep) 78 82 throw (utility::IO_error,std::exception) 79 : blas_result_(NULL), view_(NULL) 83 : blas_result_(NULL), view_(NULL), view_const_(NULL) 80 84 { 81 85 // read the data file and store in stl vectors (dynamically … … 133 137 is.clear(std::ios::goodbit); 134 138 // convert the data to a gsl matrix 135 m_ = gsl_matrix_alloc ( nof_rows, nof_columns );139 proxy_m_ = m_ = gsl_matrix_alloc ( nof_rows, nof_columns ); 136 140 if (!m_) 137 141 throw utility::GSL_error("matrix::matrix failed to allocate memory"); … … 149 153 if (view_) 150 154 delete view_; 151 else {152 if (blas_result_)153 gsl_matrix_free(blas_result_);154 if (m_)155 gsl_matrix_free(m_);156 }155 else if (view_const_) 156 delete view_const_; 157 else if (m_) 158 gsl_matrix_free(m_); 159 if (blas_result_) 160 gsl_matrix_free(blas_result_); 157 161 blas_result_=NULL; 158 162 m_=NULL; … … 160 164 161 165 166 const matrix& matrix::clone(const matrix& other) 167 { 168 if (this!=&other) { 169 170 if (view_) 171 delete view_; 172 else if (view_const_) 173 delete view_const_; 174 else if (m_) 175 gsl_matrix_free(m_); 176 if (blas_result_) { 177 gsl_matrix_free(blas_result_); 178 blas_result_=NULL; 179 } 180 proxy_m_=m_=NULL; 181 182 if (other.view_) { 183 view_ = new gsl_matrix_view(*other.view_); 184 proxy_m_ = m_ = &(view_>matrix); 185 } 186 else if (other.view_const_) { 187 view_const_ = new gsl_matrix_const_view(*other.view_const_); 188 proxy_m_ = &(view_const_>matrix); 189 } 190 else if (other.m_) 191 proxy_m_ = m_ = other.create_gsl_matrix_copy(); 192 193 } 194 return *this; 195 } 196 197 162 198 size_t matrix::columns(void) const 163 199 { 164 if (! m_)200 if (!proxy_m_) 165 201 return 0; 166 return m_>size2;202 return proxy_m_>size2; 167 203 } 168 204 … … 173 209 if (!m) 174 210 throw utility::GSL_error("matrix::create_gsl_matrix_copy failed to allocate memory"); 175 if (gsl_matrix_memcpy(m, m_))211 if (gsl_matrix_memcpy(m,proxy_m_)) 176 212 throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mismatch"); 177 213 return m; … … 179 215 180 216 181 void matrix::div_elements(const matrix& b) 182 { 183 int status=gsl_matrix_div_elements(m_, b.m_); 217 void matrix::div(const matrix& other) 218 { 219 assert(m_); 220 int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p()); 184 221 if (status) 185 222 throw utility::GSL_error(std::string("matrix::div_elements",status)); … … 205 242 const gsl_matrix* matrix::gsl_matrix_p(void) const 206 243 { 244 return proxy_m_; 245 } 246 247 248 gsl_matrix* matrix::gsl_matrix_p(void) 249 { 207 250 return m_; 208 251 } 209 252 210 253 211 gsl_matrix* matrix::gsl_matrix_p(void)212 {213 return m_;214 }215 216 217 254 bool matrix::isview(void) const 218 255 { 219 return view_; 220 } 221 222 223 void matrix::mul_elements(const matrix& b) 224 { 225 int status=gsl_matrix_mul_elements(m_, b.m_); 256 return view_  view_const_; 257 } 258 259 260 void matrix::mul(const matrix& other) 261 { 262 assert(m_); 263 int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p()); 226 264 if (status) 227 265 throw utility::GSL_error(std::string("matrix::mul_elements",status)); 228 229 266 } 230 267 … … 232 269 size_t matrix::rows(void) const 233 270 { 234 if (! m_)271 if (!proxy_m_) 235 272 return 0; 236 return m_>size1; 237 } 238 239 240 void matrix::set(const matrix& mat) 241 { 242 if (gsl_matrix_memcpy(m_, mat.m_)) 273 return proxy_m_>size1; 274 } 275 276 277 void matrix::set(const matrix& other) 278 { 279 assert(m_); 280 if (gsl_matrix_memcpy(m_, other.gsl_matrix_p())) 243 281 throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mismatch"); 244 282 } … … 247 285 void matrix::set_all(const double value) 248 286 { 287 assert(m_); 249 288 gsl_matrix_set_all(m_, value); 250 289 } … … 253 292 void matrix::set_column(const size_t column, const vector& vec) 254 293 { 294 assert(m_); 255 295 int status=gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); 256 296 if (status) … … 261 301 void matrix::set_row(const size_t row, const vector& vec) 262 302 { 303 assert(m_); 263 304 int status=gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); 264 305 if (status) … … 269 310 void matrix::swap_columns(const size_t i, const size_t j) 270 311 { 312 assert(m_); 271 313 int status=gsl_matrix_swap_columns(m_, i, j); 272 314 if (status) … … 277 319 void matrix::swap_rowcol(const size_t i, const size_t j) 278 320 { 321 assert(m_); 279 322 int status=gsl_matrix_swap_rowcol(m_, i, j); 280 323 if (status) … … 285 328 void matrix::swap_rows(const size_t i, const size_t j) 286 329 { 330 assert(m_); 287 331 int status=gsl_matrix_swap_rows(m_, i, j); 288 332 if (status) … … 293 337 void matrix::transpose(void) 294 338 { 339 assert(m_); 295 340 if (columns()==rows()) 296 341 gsl_matrix_transpose(m_); // this never fails … … 302 347 gsl_matrix_transpose_memcpy(transposed,m_); 303 348 gsl_matrix_free(m_); 304 m_=transposed;349 proxy_m_ = m_ = transposed; 305 350 if (blas_result_) { 306 351 gsl_matrix_free(blas_result_); … … 313 358 double& matrix::operator()(size_t row, size_t column) 314 359 { 360 assert(m_); 315 361 double* d=gsl_matrix_ptr(m_, row, column); 316 362 if (!d) … … 322 368 const double& matrix::operator()(size_t row, size_t column) const 323 369 { 324 const double* d=gsl_matrix_const_ptr( m_, row, column);370 const double* d=gsl_matrix_const_ptr(proxy_m_, row, column); 325 371 if (!d) 326 372 throw utility::GSL_error("matrix::operator()",GSL_EINVAL); … … 343 389 const matrix& matrix::operator=( const matrix& other ) 344 390 { 345 if ( this != &other ) { 346 if (view_) 347 delete view_; 348 else if (m_) 349 gsl_matrix_free(m_); 350 m_ = other.create_gsl_matrix_copy(); 351 } 352 return *this; 353 } 354 355 356 const matrix& matrix::operator+=(const matrix& m) 357 { 358 int status=gsl_matrix_add(m_, m.m_); 391 set(other); 392 return *this; 393 } 394 395 396 const matrix& matrix::operator+=(const matrix& other) 397 { 398 assert(m_); 399 int status=gsl_matrix_add(m_, other.proxy_m_); 359 400 if (status) 360 401 throw utility::GSL_error(std::string("matrix::operator+=", status)); … … 365 406 const matrix& matrix::operator+=(const double d) 366 407 { 408 assert(m_); 367 409 gsl_matrix_add_constant(m_, d); 368 410 return *this; … … 370 412 371 413 372 const matrix& matrix::operator=(const matrix& m) 373 { 374 int status=gsl_matrix_sub(m_, m.m_); 414 const matrix& matrix::operator=(const matrix& other) 415 { 416 assert(m_); 417 int status=gsl_matrix_sub(m_, other.proxy_m_); 375 418 if (status) 376 419 throw utility::GSL_error(std::string("matrix::operator=", status)); … … 381 424 const matrix& matrix::operator=(const double d) 382 425 { 426 assert(m_); 383 427 gsl_matrix_add_constant(m_, d); 384 428 return *this; … … 388 432 const matrix& matrix::operator*=(const matrix& other) 389 433 { 434 assert(m_); 390 435 if (!blas_result_) { 391 436 blas_result_ = gsl_matrix_alloc(rows(),other.columns()); … … 393 438 throw utility::GSL_error("matrix::operator*= failed to allocate memory"); 394 439 } 395 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other. m_, 0.0,440 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.proxy_m_, 0.0, 396 441 blas_result_); 397 442 gsl_matrix* tmp=m_; 398 m_=blas_result_;443 proxy_m_ = m_ = blas_result_; 399 444 blas_result_=tmp; 400 445 return *this; … … 404 449 const matrix& matrix::operator*=(const double d) 405 450 { 451 assert(m_); 406 452 gsl_matrix_scale(m_, d); 407 453 return *this; … … 409 455 410 456 411 bool isnull(const matrix& m)412 { 413 return gsl_matrix_isnull( m.gsl_matrix_p());414 } 415 416 417 double max(const matrix& m)418 { 419 return gsl_matrix_max( m.gsl_matrix_p());420 } 421 422 423 double min(const matrix& m)424 { 425 return gsl_matrix_min( m.gsl_matrix_p());426 } 427 428 429 void minmax_index(const matrix& m,457 bool isnull(const matrix& other) 458 { 459 return gsl_matrix_isnull(other.gsl_matrix_p()); 460 } 461 462 463 double max(const matrix& other) 464 { 465 return gsl_matrix_max(other.gsl_matrix_p()); 466 } 467 468 469 double min(const matrix& other) 470 { 471 return gsl_matrix_min(other.gsl_matrix_p()); 472 } 473 474 475 void minmax_index(const matrix& other, 430 476 std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max) 431 477 { 432 gsl_matrix_minmax_index( m.gsl_matrix_p(), &min.first, &min.second,478 gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second, 433 479 &max.first, &max.second); 434 480 } … … 440 486 size_t columns=templat.columns(); 441 487 if (rows!=flag.rows() && columns!=flag.columns()) 442 flag =matrix(rows,columns,1.0);488 flag.clone(matrix(rows,columns,1.0)); 443 489 else 444 490 flag.set_all(1.0); … … 456 502 void swap(matrix& a, matrix& b) 457 503 { 504 assert(v.gsl_matrix_p()); assert(w.gsl_matrix_p()); 458 505 int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p()); 459 506 if (status) 
trunk/yat/utility/matrix.h
r788 r792 64 64 /// supported are no binary file support and views on arrays, 65 65 /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and 66 /// superdiagonals. If there is an interest from the user community, 67 /// the omitted functionality could be included. 66 /// superdiagonals. 68 67 /// 69 68 class matrix … … 143 142 ~matrix(void); 144 143 144 /** 145 \brief Make a copy of \a other. 146 147 This function will make a deep copy of \a other. Memory is 148 resized and view state is changed if needed. 149 */ 150 const matrix& clone(const matrix& other); 151 145 152 /// 146 153 /// @return The number of columns in the matrix. … … 155 162 \throw GSL_error if dimensions mismatch. 156 163 */ 157 void div _elements(const matrix& b);164 void div(const matrix& b); 158 165 159 166 /** … … 193 200 \throw GSL_error if dimensions mismatch. 194 201 */ 195 void mul _elements(const matrix& b);202 void mul(const matrix& b); 196 203 197 204 /// … … 320 327 \brief The assignment operator. 321 328 322 There is no requirements on dimensions, i.e. the matrix is 323 remapped in memory if necessary. This implies that in general 324 views cannot be assigned using this operator. Views will be 325 mutated into normal matrices. The only exception to this 326 behaviour on views is when selfassignemnt is done, since 327 selfassignment is ignored. 329 Dimensions of the matrices must match. If the LHS matrix is a 330 view, the underlying data will be changed. 328 331 329 332 \return A const reference to the resulting matrix. … … 331 334 \see void set(const matrix&). 332 335 333 \throw A GSL_error is indirectly thrown if memory allocation 334 fails, or if dimensions mismatch. 336 \throw GSL_error if dimensions mismatch. 335 337 */ 336 338 const matrix& operator=(const matrix& other); … … 419 421 gsl_matrix* m_; 420 422 gsl_matrix_view* view_; 423 const gsl_matrix_const_view* view_const_; 424 // proxy_m_ is used to access the proper underlying gsl_matrix in 425 // all const member functions. It is not used by design for 426 // nonconst vector functions and operators. This is to make sure 427 // that runtime errors occur if a const matrix is used in an 428 // inappropriate manner such as on left hand side in assignment 429 // (remember, m_ is null for const matrix views). 430 const gsl_matrix* proxy_m_; 421 431 }; 422 432
Note: See TracChangeset
for help on using the changeset viewer.