Changeset 762
 Timestamp:
 Feb 20, 2007, 8:18:48 PM (15 years ago)
 Location:
 trunk
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test/matrix_test.cc
r680 r762 58 58 *error << "Testing matrix class" << std::endl; 59 59 bool ok = true; 60 utility::matrix unit3x3(3,3); 61 for (size_t i=0; i<unit3x3.rows(); ++i) 62 unit3x3(i,i)=1; 60 63 61 64 *error << "\tcopy constructor and operator!=" << std::endl; … … 78 81 std::remove("data/tmp_test_matrix.txt"); 79 82 80 *error << "\toperator* (double)" << std::endl;83 *error << "\toperator*=(double)" << std::endl; 81 84 utility::matrix m4(3,3,1); 82 85 m4 *= 9; … … 233 236 delete m_nan; 234 237 238 *error << "\toperator*=(matrix&)" << std::endl; 239 utility::matrix m6(unit3x3); 240 m6 *= m; 241 if (m6!=m) { 242 ok=false; 243 *error << "error operator*=(matrix) 1" << std::endl; 244 } 245 m6 *= unit3x3; 246 if (m6!=m) { 247 ok=false; 248 *error << "error operator*=(matrix) 2" << std::endl; 249 } 250 235 251 if (error!=&std::cerr) 236 252 delete error; 
trunk/yat/utility/matrix.cc
r759 r762 41 41 42 42 matrix::matrix(void) 43 : m_(NULL), view_(NULL)43 : blas_result_(NULL), m_(NULL), view_(NULL) 44 44 { 45 45 } … … 47 47 48 48 matrix::matrix(const size_t& r, const size_t& c, double init_value) 49 : m_(gsl_matrix_alloc(r,c)), view_(NULL)49 : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL) 50 50 { 51 51 if (!m_) … … 56 56 57 57 matrix::matrix(const matrix& o) 58 : m_(o.create_gsl_matrix_copy()), view_(NULL)58 : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL) 59 59 { 60 60 } … … 63 63 matrix::matrix(matrix& m, size_t offset_row, size_t offset_column, 64 64 size_t n_row, size_t n_column) 65 : blas_result_(NULL) 65 66 { 66 67 view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_, … … 76 77 matrix::matrix(std::istream& is, char sep) 77 78 throw (utility::IO_error,std::exception) 78 : view_(NULL)79 : blas_result_(NULL), view_(NULL) 79 80 { 80 81 // read the data file and store in stl vectors (dynamically … … 148 149 if (view_) 149 150 delete view_; 150 else if (m_) 151 gsl_matrix_free(m_); 151 else { 152 if (blas_result_) 153 gsl_matrix_free(blas_result_); 154 if (m_) 155 gsl_matrix_free(m_); 156 } 157 blas_result_=NULL; 152 158 m_=NULL; 153 159 } … … 371 377 gsl_matrix_free(m_); 372 378 m_=transposed; 379 if (blas_result_) { 380 gsl_matrix_free(blas_result_); 381 blas_result_=NULL; 382 } 373 383 } 374 384 } … … 441 451 const matrix& matrix::operator*=(const matrix& other) 442 452 { 443 gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns()); 444 if (!result) 445 throw utility::GSL_error("matrix::operator*= failed to allocate memory"); 446 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result); 447 gsl_matrix_free(m_); 448 m_=result; 453 if (!blas_result_) { 454 blas_result_ = gsl_matrix_alloc(rows(),other.columns()); 455 if (!blas_result_) 456 throw utility::GSL_error("matrix::operator*= failed to allocate memory"); 457 } 458 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, 459 blas_result_); 460 gsl_matrix* tmp=m_; 461 m_=blas_result_; 462 blas_result_=tmp; 449 463 return *this; 450 464 } 
trunk/yat/utility/matrix.h
r759 r762 67 67 /// the omitted functionality could be included. 68 68 /// 69 /// @todo Maybe it would be smart to create temporary objects need70 /// for BLAS in constructors. As it is now temporary objects are71 /// called before BLAS functionality iss used, cf. const matrix&72 /// matrix::operator*=(const matrix& other)73 ///74 69 class matrix 75 70 { … … 323 318 \brief Swap row \a i and column \a j. 324 319 320 The matrix must be square. 321 325 322 \throw GSL_error if either index is out of bounds, or if matrix 326 323 is not square. … … 457 454 gsl_matrix* create_gsl_matrix_copy(void) const; 458 455 456 // blas_result_ is used to temporarily store result in BLAS calls 457 // and when not NULL it should always have the same size as 458 // m_. Memory is not allocated for blas_result_ until it is 459 // needed. 460 gsl_matrix* blas_result_; 459 461 gsl_matrix* m_; 460 462 gsl_matrix_view* view_;
Note: See TracChangeset
for help on using the changeset viewer.