Changeset 1380 for trunk/yat/utility/MatrixWeighted.cc
 Timestamp:
 Jul 17, 2008, 12:36:23 AM (14 years ago)
 File:

 1 copied
Legend:
 Unmodified
 Added
 Removed

trunk/yat/utility/MatrixWeighted.cc
r1360 r1380 25 25 */ 26 26 27 #include "Matrix.h" 28 #include "Vector.h" 29 #include "VectorBase.h" 30 #include "VectorConstView.h" 31 #include "VectorView.h" 32 #include "utility.h" 27 #include "MatrixWeighted.h" 33 28 34 29 #include <cassert> 35 #include <cmath>36 #include <sstream>37 30 #include <vector> 38 39 #include <gsl/gsl_blas.h>40 31 41 32 namespace theplu { … … 44 35 45 36 46 Matrix ::Matrix(void)47 : blas_result_(NULL), m_(NULL)37 MatrixWeighted::MatrixWeighted(void) 38 : columns_(0) 48 39 { 49 40 } 50 41 51 42 52 Matrix::Matrix(const size_t& r, const size_t& c, double init_value)53 : blas_result_(NULL), m_(gsl_matrix_alloc(r,c))54 {55 if (!m_)56 throw utility::GSL_error("Matrix::Matrix failed to allocate memory");57 all(init_value);58 }59 60 61 Matrix::Matrix(const Matrix& o)62 : blas_result_(NULL), m_(o.create_gsl_matrix_copy())63 {64 }65 66 67 // Constructor that gets data from istream68 Matrix::Matrix(std::istream& is, char sep)69 throw (utility::IO_error,std::exception)70 : blas_result_(NULL)71 {72 if (!is.good())73 throw utility::IO_error("Matrix: istream is not good");74 75 // read the data file and store in stl vectors (dynamically76 // expandable)77 std::vector<std::vector<double> > data_matrix;78 unsigned int nof_columns=0;79 unsigned int nof_rows = 0;80 std::string line;81 while(getline(is, line, '\n')){82 // Ignoring empty lines83 if (!line.size()) {84 continue;85 }86 nof_rows++;87 std::vector<double> v;88 std::string element;89 std::stringstream ss(line);90 91 bool ok=true;92 while(ok) {93 if(sep=='\0')94 ok=(ss>>element);95 else96 ok=getline(ss, element, sep);97 if(!ok)98 break;99 100 if (!element.size())101 v.push_back(std::numeric_limits<double>::quiet_NaN());102 else {103 try {104 v.push_back(convert<double>(element));105 }106 catch (std::runtime_error& e) {107 std::stringstream ss(e.what());108 ss << "\nMatrix.cc: " << element109 << " is not accepted as a matrix element\n";110 throw IO_error(ss.str());111 }112 }113 }114 if(sep!='\0' && line[line.size()1]==sep) // add NaN for final separator115 v.push_back(std::numeric_limits<double>::quiet_NaN());116 if (!nof_columns)117 nof_columns=v.size();118 else if (v.size()!=nof_columns) {119 std::ostringstream s;120 s << "Matrix::Matrix(std::istream&, char) data file error: "121 << "line " << nof_rows << " has " << v.size()122 << " columns; expected " << nof_columns << " columns.";123 throw utility::IO_error(s.str());124 }125 data_matrix.push_back(v);126 }127 128 // manipulate the state of the stream to be good129 is.clear(std::ios::goodbit);130 131 // if stream was empty, create nothing132 if (!nof_columns  !nof_rows)133 return;134 135 // convert the data to a gsl matrix136 m_ = gsl_matrix_alloc ( nof_rows, nof_columns );137 if (!m_)138 throw utility::GSL_error("Matrix::Matrix failed to allocate memory");139 140 // if gsl error handler disabled, out of bounds index will not141 // abort the program.142 for(size_t i=0;i<nof_rows;i++)143 for(size_t j=0;j<nof_columns;j++)144 gsl_matrix_set( m_, i, j, data_matrix[i][j] );145 }146 147 148 Matrix::~Matrix(void)149 {150 delete_allocated_memory();151 if (blas_result_)152 gsl_matrix_free(blas_result_);153 blas_result_=NULL;154 }155 156 157 void Matrix::all(const double value)158 {159 assert(m_);160 gsl_matrix_set_all(m_, value);161 }162 163 164 Matrix::iterator Matrix::begin(void)165 {166 return iterator(&(*this)(0,0), 1);167 }168 169 170 Matrix::const_iterator Matrix::begin(void) const171 {172 return const_iterator(&(*this)(0,0), 1);173 }174 175 176 Matrix::column_iterator Matrix::begin_column(size_t i)177 {178 return iterator(&(*this)(0,i), this>columns());179 }180 181 182 Matrix::const_column_iterator Matrix::begin_column(size_t i) const183 {184 return const_iterator(&(*this)(0,i), this>columns());185 }186 187 188 Matrix::row_iterator Matrix::begin_row(size_t i)189 {190 return iterator(&(*this)(i,0), 1);191 }192 193 194 Matrix::const_row_iterator Matrix::begin_row(size_t i) const195 {196 return const_iterator(&(*this)(i,0), 1);197 }198 199 200 VectorView Matrix::column_view(size_t col)201 {202 VectorView res(*this, col, false);203 return res;204 }205 206 207 const VectorConstView Matrix::column_const_view(size_t col) const208 {209 return VectorConstView(*this, col, false);210 }211 212 213 size_t Matrix::columns(void) const214 {215 return (m_ ? m_>size2 : 0);216 }217 218 219 gsl_matrix* Matrix::create_gsl_matrix_copy(void) const220 {221 gsl_matrix* m = gsl_matrix_alloc(rows(),columns());222 if (!m)223 throw utility::GSL_error("Matrix::create_gsl_matrix_copy failed to allocate memory");224 if (gsl_matrix_memcpy(m,m_))225 throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mismatch");226 return m;227 }228 229 230 void Matrix::delete_allocated_memory(void)231 {232 if (m_)233 gsl_matrix_free(m_);234 m_=NULL;235 }236 237 238 void Matrix::div(const Matrix& other)239 {240 assert(m_);241 int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());242 if (status)243 throw utility::GSL_error(std::string("matrix::div_elements",status));244 }245 246 247 Matrix::iterator Matrix::end(void)248 {249 return iterator(&(*this)(0,0)+rows()*columns(), 1);250 }251 252 253 Matrix::const_iterator Matrix::end(void) const254 {255 return const_iterator(&(*this)(0,0)+rows()*columns(), 1);256 }257 258 259 Matrix::column_iterator Matrix::end_column(size_t i)260 {261 return column_iterator(&(*this)(0,i)+rows()*columns(), this>columns());262 }263 264 265 Matrix::const_column_iterator Matrix::end_column(size_t i) const266 {267 return const_column_iterator(&(*this)(0,i)+rows()*columns(),this>columns());268 }269 270 271 Matrix::row_iterator Matrix::end_row(size_t i)272 {273 return row_iterator(&(*this)(i,0)+columns(), 1);274 }275 276 277 Matrix::const_row_iterator Matrix::end_row(size_t i) const278 {279 return const_row_iterator(&(*this)(i,0)+columns(), 1);280 }281 282 283 bool Matrix::equal(const Matrix& other, const double d) const284 {285 if (this==&other)286 return true;287 if (columns()!=other.columns()  rows()!=other.rows())288 return false;289 for (size_t i=0; i<rows(); i++)290 for (size_t j=0; j<columns(); j++)291 // The two last condition checks are needed for NaN detection292 if (fabs( (*this)(i,j)other(i,j) ) > d 293 (*this)(i,j)!=(*this)(i,j)  other(i,j)!=other(i,j))294 return false;295 return true;296 }297 298 299 const gsl_matrix* Matrix::gsl_matrix_p(void) const300 {301 return m_;302 }303 304 305 gsl_matrix* Matrix::gsl_matrix_p(void)306 {307 return m_;308 }309 310 311 void Matrix::mul(const Matrix& other)312 {313 assert(m_);314 int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());315 if (status)316 throw utility::GSL_error(std::string("Matrix::mul_elements",status));317 }318 319 320 void Matrix::resize(size_t r, size_t c, double init_value)321 {322 delete_allocated_memory();323 324 assert( (r&&c)  (!r&&!c) );325 if (r && c){326 m_ = gsl_matrix_alloc(r,c);327 if (!m_)328 throw utility::GSL_error("Matrix::Matrix failed to allocate memory");329 all(init_value);330 }331 332 // no need to delete blas_result_ if the number of rows fit, it333 // may be useful later.334 if (blas_result_ && (blas_result_>size1!=rows())) {335 gsl_matrix_free(blas_result_);336 blas_result_=NULL;337 }338 }339 340 341 size_t Matrix::rows(void) const342 {343 return (m_ ? m_>size1 : 0);344 }345 346 347 const VectorConstView Matrix::row_const_view(size_t col) const348 {349 return VectorConstView(*this, col, true);350 }351 352 353 VectorView Matrix::row_view(size_t row)354 {355 VectorView res(*this, row, true);356 return res;357 }358 359 360 void Matrix::swap_columns(const size_t i, const size_t j)361 {362 assert(m_);363 int status=gsl_matrix_swap_columns(m_, i, j);364 if (status)365 throw utility::GSL_error(std::string("Matrix::swap_columns",status));366 }367 368 369 void Matrix::swap_rowcol(const size_t i, const size_t j)370 {371 assert(m_);372 int status=gsl_matrix_swap_rowcol(m_, i, j);373 if (status)374 throw utility::GSL_error(std::string("Matrix::swap_rowcol",status));375 }376 377 378 void Matrix::swap_rows(const size_t i, const size_t j)379 {380 assert(m_);381 int status=gsl_matrix_swap_rows(m_, i, j);382 if (status)383 throw utility::GSL_error(std::string("Matrix::swap_rows",status));384 }385 386 387 void Matrix::transpose(void)388 {389 assert(m_);390 if (columns()==rows())391 gsl_matrix_transpose(m_); // this never fails392 else {393 gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());394 if (!transposed)395 throw utility::GSL_error("Matrix::transpose failed to allocate memory");396 // next line never fails if allocation above succeeded.397 gsl_matrix_transpose_memcpy(transposed,m_);398 gsl_matrix_free(m_);399 m_ = transposed;400 if (blas_result_) {401 gsl_matrix_free(blas_result_);402 blas_result_=NULL;403 }404 }405 }406 407 408 double& Matrix::operator()(size_t row, size_t column)409 {410 assert(m_);411 assert(row<rows());412 assert(column<columns());413 double* d=gsl_matrix_ptr(m_, row, column);414 if (!d)415 throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);416 return *d;417 }418 419 420 const double& Matrix::operator()(size_t row, size_t column) const421 {422 assert(row<rows());423 assert(column<columns());424 const double* d=gsl_matrix_const_ptr(m_, row, column);425 if (!d)426 throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);427 return *d;428 }429 430 431 bool Matrix::operator==(const Matrix& other) const432 {433 return equal(other);434 }435 436 437 bool Matrix::operator!=(const Matrix& other) const438 {439 return !equal(other);440 }441 442 443 const Matrix& Matrix::operator=( const Matrix& other )444 {445 assert(other.m_);446 if (this!=&other) {447 if ( !m_  (other.m_>size1!=m_>size1)  (other.m_>size2!=m_>size2) )448 resize(other.m_>size1,other.m_>size2);449 if (gsl_matrix_memcpy(m_, other.gsl_matrix_p()))450 throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mismatch");451 }452 return *this;453 }454 455 456 const Matrix& Matrix::operator+=(const Matrix& other)457 {458 assert(m_);459 int status=gsl_matrix_add(m_, other.m_);460 if (status)461 throw utility::GSL_error(std::string("Matrix::operator+=", status));462 return *this;463 }464 465 466 const Matrix& Matrix::operator+=(const double d)467 {468 assert(m_);469 gsl_matrix_add_constant(m_, d);470 return *this;471 }472 473 474 const Matrix& Matrix::operator=(const Matrix& other)475 {476 assert(m_);477 int status=gsl_matrix_sub(m_, other.m_);478 if (status)479 throw utility::GSL_error(std::string("Matrix::operator=", status));480 return *this;481 }482 483 484 const Matrix& Matrix::operator=(const double d)485 {486 assert(m_);487 gsl_matrix_add_constant(m_, d);488 return *this;489 }490 491 492 const Matrix& Matrix::operator*=(const Matrix& other)493 {494 assert(m_);495 if ( blas_result_ && ((blas_result_>size1!=rows()) 496 (blas_result_>size2!=other.columns())) ) {497 gsl_matrix_free(blas_result_);498 blas_result_=NULL;499 }500 if (!blas_result_) {501 blas_result_ = gsl_matrix_alloc(rows(),other.columns());502 if (!blas_result_)503 throw utility::GSL_error("Matrix::operator*= failed to allocate memory");504 }505 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, blas_result_);506 gsl_matrix* tmp=m_;507 m_ = blas_result_;508 blas_result_=tmp;509 return *this;510 }511 512 513 const Matrix& Matrix::operator*=(const double d)514 {515 assert(m_);516 gsl_matrix_scale(m_, d);517 return *this;518 }519 520 521 bool isnull(const Matrix& other)522 {523 return gsl_matrix_isnull(other.gsl_matrix_p());524 }525 526 527 double max(const Matrix& other)528 {529 return gsl_matrix_max(other.gsl_matrix_p());530 }531 532 533 double min(const Matrix& other)534 {535 return gsl_matrix_min(other.gsl_matrix_p());536 }537 538 539 void minmax_index(const Matrix& other,540 std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)541 {542 gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,543 &max.first, &max.second);544 }545 546 547 bool nan(const Matrix& templat, Matrix& flag)548 {549 size_t rows=templat.rows();550 size_t columns=templat.columns();551 if (rows!=flag.rows() && columns!=flag.columns())552 flag.resize(rows,columns,1.0);553 else554 flag.all(1.0);555 bool nan=false;556 for (size_t i=0; i<rows; i++)557 for (size_t j=0; j<columns; j++)558 if (std::isnan(templat(i,j))) {559 flag(i,j)=0;560 nan=true;561 }562 return nan;563 }564 565 566 void swap(Matrix& a, Matrix& b)567 {568 assert(a.gsl_matrix_p()); assert(b.gsl_matrix_p());569 int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());570 if (status)571 throw utility::GSL_error(std::string("swap(Matrix&,Matrix&)",status));572 }573 574 575 std::ostream& operator<<(std::ostream& s, const Matrix& m)576 {577 s.setf(std::ios::dec);578 s.precision(12);579 for(size_t i=0, j=0; i<m.rows(); i++)580 for (j=0; j<m.columns(); j++) {581 s << m(i,j);582 if (j<m.columns()1)583 s << s.fill();584 else if (i<m.rows()1)585 s << "\n";586 }587 return s;588 }589 590 591 Vector operator*(const Matrix& m, const VectorBase& v)592 {593 utility::Vector res(m.rows());594 for (size_t i=0; i<res.size(); ++i)595 res(i) = VectorConstView(m,i) * v;596 return res;597 }598 599 600 Vector operator*(const VectorBase& v, const Matrix& m)601 {602 utility::Vector res(m.columns());603 for (size_t i=0; i<res.size(); ++i)604 res(i) = v * VectorConstView(m,i,false);605 return res;606 }607 608 43 }}} // of namespace utility, yat and thep
Note: See TracChangeset
for help on using the changeset viewer.