[42] | 1 | // $Id: Matrix.cc 1598 2008-10-24 21:12:37Z peter $ |
---|
[12] | 2 | |
---|
[570] | 3 | /* |
---|
| 4 | Copyright (C) 2003 Daniel Dalevi, Peter Johansson |
---|
[831] | 5 | Copyright (C) 2004 Jari Häkkinen, Peter Johansson |
---|
[1275] | 6 | Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér |
---|
[1260] | 7 | Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson |
---|
[570] | 8 | |
---|
[1437] | 9 | This file is part of the yat library, http://dev.thep.lu.se/yat |
---|
[570] | 10 | |
---|
[675] | 11 | The yat library is free software; you can redistribute it and/or |
---|
| 12 | modify it under the terms of the GNU General Public License as |
---|
[1486] | 13 | published by the Free Software Foundation; either version 3 of the |
---|
[675] | 14 | License, or (at your option) any later version. |
---|
[570] | 15 | |
---|
[675] | 16 | The yat library is distributed in the hope that it will be useful, |
---|
| 17 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 18 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
---|
[570] | 19 | General Public License for more details. |
---|
| 20 | |
---|
| 21 | You should have received a copy of the GNU General Public License |
---|
[1487] | 22 | along with yat. If not, see <http://www.gnu.org/licenses/>. |
---|
[570] | 23 | */ |
---|
| 24 | |
---|
[1121] | 25 | #include "Matrix.h" |
---|
| 26 | #include "Vector.h" |
---|
[1017] | 27 | #include "VectorBase.h" |
---|
[1027] | 28 | #include "VectorConstView.h" |
---|
[1015] | 29 | #include "VectorView.h" |
---|
[680] | 30 | #include "utility.h" |
---|
[273] | 31 | |
---|
[792] | 32 | #include <cassert> |
---|
[293] | 33 | #include <cmath> |
---|
[42] | 34 | #include <sstream> |
---|
[17] | 35 | #include <vector> |
---|
| 36 | |
---|
[42] | 37 | #include <gsl/gsl_blas.h> |
---|
| 38 | |
---|
| 39 | namespace theplu { |
---|
[680] | 40 | namespace yat { |
---|
[616] | 41 | namespace utility { |
---|
[12] | 42 | |
---|
| 43 | |
---|
[1121] | 44 | Matrix::Matrix(void) |
---|
[1098] | 45 | : blas_result_(NULL), m_(NULL) |
---|
[703] | 46 | { |
---|
| 47 | } |
---|
[12] | 48 | |
---|
[703] | 49 | |
---|
[1121] | 50 | Matrix::Matrix(const size_t& r, const size_t& c, double init_value) |
---|
[1598] | 51 | : blas_result_(NULL), m_(NULL) |
---|
[703] | 52 | { |
---|
[1598] | 53 | resize(r,c,init_value); |
---|
[703] | 54 | } |
---|
| 55 | |
---|
| 56 | |
---|
[1121] | 57 | Matrix::Matrix(const Matrix& o) |
---|
[1098] | 58 | : blas_result_(NULL), m_(o.create_gsl_matrix_copy()) |
---|
[703] | 59 | { |
---|
| 60 | } |
---|
| 61 | |
---|
| 62 | |
---|
[42] | 63 | // Constructor that gets data from istream |
---|
[1121] | 64 | Matrix::Matrix(std::istream& is, char sep) |
---|
[434] | 65 | throw (utility::IO_error,std::exception) |
---|
[1098] | 66 | : blas_result_(NULL) |
---|
[42] | 67 | { |
---|
[1147] | 68 | if (!is.good()) |
---|
| 69 | throw utility::IO_error("Matrix: istream is not good"); |
---|
| 70 | |
---|
[42] | 71 | // read the data file and store in stl vectors (dynamically |
---|
| 72 | // expandable) |
---|
| 73 | std::vector<std::vector<double> > data_matrix; |
---|
[1271] | 74 | unsigned int nof_columns=0; |
---|
| 75 | unsigned int nof_rows = 0; |
---|
[434] | 76 | std::string line; |
---|
| 77 | while(getline(is, line, '\n')){ |
---|
[128] | 78 | // Ignoring empty lines |
---|
[434] | 79 | if (!line.size()) { |
---|
[128] | 80 | continue; |
---|
| 81 | } |
---|
[434] | 82 | nof_rows++; |
---|
| 83 | std::vector<double> v; |
---|
| 84 | std::string element; |
---|
| 85 | std::stringstream ss(line); |
---|
[439] | 86 | |
---|
[434] | 87 | bool ok=true; |
---|
| 88 | while(ok) { |
---|
| 89 | if(sep=='\0') |
---|
| 90 | ok=(ss>>element); |
---|
| 91 | else |
---|
| 92 | ok=getline(ss, element, sep); |
---|
| 93 | if(!ok) |
---|
| 94 | break; |
---|
[567] | 95 | |
---|
[1313] | 96 | if (!element.size()) |
---|
[434] | 97 | v.push_back(std::numeric_limits<double>::quiet_NaN()); |
---|
| 98 | else { |
---|
[1313] | 99 | try { |
---|
| 100 | v.push_back(convert<double>(element)); |
---|
| 101 | } |
---|
| 102 | catch (std::runtime_error& e) { |
---|
| 103 | std::stringstream ss(e.what()); |
---|
| 104 | ss << "\nMatrix.cc: " << element |
---|
| 105 | << " is not accepted as a matrix element\n"; |
---|
| 106 | throw IO_error(ss.str()); |
---|
| 107 | } |
---|
[434] | 108 | } |
---|
| 109 | } |
---|
[439] | 110 | if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator |
---|
| 111 | v.push_back(std::numeric_limits<double>::quiet_NaN()); |
---|
[420] | 112 | if (!nof_columns) |
---|
[42] | 113 | nof_columns=v.size(); |
---|
[420] | 114 | else if (v.size()!=nof_columns) { |
---|
| 115 | std::ostringstream s; |
---|
[1121] | 116 | s << "Matrix::Matrix(std::istream&, char) data file error: " |
---|
[444] | 117 | << "line " << nof_rows << " has " << v.size() |
---|
[420] | 118 | << " columns; expected " << nof_columns << " columns."; |
---|
| 119 | throw utility::IO_error(s.str()); |
---|
[42] | 120 | } |
---|
| 121 | data_matrix.push_back(v); |
---|
| 122 | } |
---|
[273] | 123 | |
---|
[42] | 124 | // manipulate the state of the stream to be good |
---|
| 125 | is.clear(std::ios::goodbit); |
---|
[1147] | 126 | |
---|
| 127 | // if stream was empty, create nothing |
---|
| 128 | if (!nof_columns || !nof_rows) |
---|
| 129 | return; |
---|
| 130 | |
---|
[42] | 131 | // convert the data to a gsl matrix |
---|
[1098] | 132 | m_ = gsl_matrix_alloc ( nof_rows, nof_columns ); |
---|
[753] | 133 | if (!m_) |
---|
[1121] | 134 | throw utility::GSL_error("Matrix::Matrix failed to allocate memory"); |
---|
[754] | 135 | |
---|
| 136 | // if gsl error handler disabled, out of bounds index will not |
---|
| 137 | // abort the program. |
---|
[1271] | 138 | for(size_t i=0;i<nof_rows;i++) |
---|
| 139 | for(size_t j=0;j<nof_columns;j++) |
---|
[42] | 140 | gsl_matrix_set( m_, i, j, data_matrix[i][j] ); |
---|
| 141 | } |
---|
[12] | 142 | |
---|
| 143 | |
---|
[1121] | 144 | Matrix::~Matrix(void) |
---|
[42] | 145 | { |
---|
[810] | 146 | delete_allocated_memory(); |
---|
[792] | 147 | if (blas_result_) |
---|
| 148 | gsl_matrix_free(blas_result_); |
---|
[1098] | 149 | blas_result_=NULL; |
---|
[42] | 150 | } |
---|
[12] | 151 | |
---|
| 152 | |
---|
[1121] | 153 | void Matrix::all(const double value) |
---|
[1098] | 154 | { |
---|
| 155 | assert(m_); |
---|
| 156 | gsl_matrix_set_all(m_, value); |
---|
| 157 | } |
---|
| 158 | |
---|
| 159 | |
---|
[1121] | 160 | Matrix::iterator Matrix::begin(void) |
---|
[1064] | 161 | { |
---|
| 162 | return iterator(&(*this)(0,0), 1); |
---|
| 163 | } |
---|
| 164 | |
---|
| 165 | |
---|
[1121] | 166 | Matrix::const_iterator Matrix::begin(void) const |
---|
[1064] | 167 | { |
---|
| 168 | return const_iterator(&(*this)(0,0), 1); |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | |
---|
[1121] | 172 | Matrix::column_iterator Matrix::begin_column(size_t i) |
---|
[1064] | 173 | { |
---|
| 174 | return iterator(&(*this)(0,i), this->columns()); |
---|
| 175 | } |
---|
| 176 | |
---|
| 177 | |
---|
[1121] | 178 | Matrix::const_column_iterator Matrix::begin_column(size_t i) const |
---|
[1064] | 179 | { |
---|
| 180 | return const_iterator(&(*this)(0,i), this->columns()); |
---|
| 181 | } |
---|
| 182 | |
---|
| 183 | |
---|
[1121] | 184 | Matrix::row_iterator Matrix::begin_row(size_t i) |
---|
[1064] | 185 | { |
---|
| 186 | return iterator(&(*this)(i,0), 1); |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | |
---|
[1121] | 190 | Matrix::const_row_iterator Matrix::begin_row(size_t i) const |
---|
[1064] | 191 | { |
---|
| 192 | return const_iterator(&(*this)(i,0), 1); |
---|
| 193 | } |
---|
| 194 | |
---|
| 195 | |
---|
[1121] | 196 | VectorView Matrix::column_view(size_t col) |
---|
[792] | 197 | { |
---|
[1098] | 198 | VectorView res(*this, col, false); |
---|
| 199 | return res; |
---|
| 200 | } |
---|
[792] | 201 | |
---|
| 202 | |
---|
[1121] | 203 | const VectorConstView Matrix::column_const_view(size_t col) const |
---|
[1098] | 204 | { |
---|
| 205 | return VectorConstView(*this, col, false); |
---|
| 206 | } |
---|
[792] | 207 | |
---|
| 208 | |
---|
[1121] | 209 | size_t Matrix::columns(void) const |
---|
[716] | 210 | { |
---|
[1098] | 211 | return (m_ ? m_->size2 : 0); |
---|
[716] | 212 | } |
---|
| 213 | |
---|
| 214 | |
---|
[1121] | 215 | gsl_matrix* Matrix::create_gsl_matrix_copy(void) const |
---|
[716] | 216 | { |
---|
[1598] | 217 | if (!m_) |
---|
| 218 | return NULL; |
---|
[716] | 219 | gsl_matrix* m = gsl_matrix_alloc(rows(),columns()); |
---|
[753] | 220 | if (!m) |
---|
[1121] | 221 | throw utility::GSL_error("Matrix::create_gsl_matrix_copy failed to allocate memory"); |
---|
[1098] | 222 | if (gsl_matrix_memcpy(m,m_)) |
---|
[1121] | 223 | throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match"); |
---|
[716] | 224 | return m; |
---|
| 225 | } |
---|
| 226 | |
---|
| 227 | |
---|
[1121] | 228 | void Matrix::delete_allocated_memory(void) |
---|
[810] | 229 | { |
---|
[1098] | 230 | if (m_) |
---|
[810] | 231 | gsl_matrix_free(m_); |
---|
[1098] | 232 | m_=NULL; |
---|
[810] | 233 | } |
---|
| 234 | |
---|
| 235 | |
---|
[1121] | 236 | void Matrix::div(const Matrix& other) |
---|
[716] | 237 | { |
---|
[792] | 238 | assert(m_); |
---|
| 239 | int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p()); |
---|
[755] | 240 | if (status) |
---|
| 241 | throw utility::GSL_error(std::string("matrix::div_elements",status)); |
---|
[716] | 242 | } |
---|
| 243 | |
---|
| 244 | |
---|
[1121] | 245 | Matrix::iterator Matrix::end(void) |
---|
[1064] | 246 | { |
---|
| 247 | return iterator(&(*this)(0,0)+rows()*columns(), 1); |
---|
| 248 | } |
---|
| 249 | |
---|
| 250 | |
---|
[1121] | 251 | Matrix::const_iterator Matrix::end(void) const |
---|
[1064] | 252 | { |
---|
| 253 | return const_iterator(&(*this)(0,0)+rows()*columns(), 1); |
---|
| 254 | } |
---|
| 255 | |
---|
| 256 | |
---|
[1121] | 257 | Matrix::column_iterator Matrix::end_column(size_t i) |
---|
[1064] | 258 | { |
---|
[1103] | 259 | return column_iterator(&(*this)(0,i)+rows()*columns(), this->columns()); |
---|
[1064] | 260 | } |
---|
| 261 | |
---|
| 262 | |
---|
[1121] | 263 | Matrix::const_column_iterator Matrix::end_column(size_t i) const |
---|
[1064] | 264 | { |
---|
[1103] | 265 | return const_column_iterator(&(*this)(0,i)+rows()*columns(),this->columns()); |
---|
[1064] | 266 | } |
---|
| 267 | |
---|
| 268 | |
---|
[1121] | 269 | Matrix::row_iterator Matrix::end_row(size_t i) |
---|
[1064] | 270 | { |
---|
[1103] | 271 | return row_iterator(&(*this)(i,0)+columns(), 1); |
---|
[1064] | 272 | } |
---|
| 273 | |
---|
| 274 | |
---|
[1121] | 275 | Matrix::const_row_iterator Matrix::end_row(size_t i) const |
---|
[1064] | 276 | { |
---|
[1103] | 277 | return const_row_iterator(&(*this)(i,0)+columns(), 1); |
---|
[1064] | 278 | } |
---|
| 279 | |
---|
| 280 | |
---|
[1121] | 281 | bool Matrix::equal(const Matrix& other, const double d) const |
---|
[272] | 282 | { |
---|
[788] | 283 | if (this==&other) |
---|
| 284 | return true; |
---|
[272] | 285 | if (columns()!=other.columns() || rows()!=other.rows()) |
---|
| 286 | return false; |
---|
[420] | 287 | for (size_t i=0; i<rows(); i++) |
---|
| 288 | for (size_t j=0; j<columns(); j++) |
---|
| 289 | // The two last condition checks are needed for NaN detection |
---|
| 290 | if (fabs( (*this)(i,j)-other(i,j) ) > d || |
---|
[293] | 291 | (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j)) |
---|
[272] | 292 | return false; |
---|
| 293 | return true; |
---|
| 294 | } |
---|
[12] | 295 | |
---|
| 296 | |
---|
[1121] | 297 | const gsl_matrix* Matrix::gsl_matrix_p(void) const |
---|
[716] | 298 | { |
---|
[1098] | 299 | return m_; |
---|
[716] | 300 | } |
---|
[12] | 301 | |
---|
[716] | 302 | |
---|
[1121] | 303 | gsl_matrix* Matrix::gsl_matrix_p(void) |
---|
[42] | 304 | { |
---|
[716] | 305 | return m_; |
---|
[42] | 306 | } |
---|
[12] | 307 | |
---|
| 308 | |
---|
[1121] | 309 | void Matrix::mul(const Matrix& other) |
---|
[716] | 310 | { |
---|
[792] | 311 | assert(m_); |
---|
| 312 | int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p()); |
---|
[755] | 313 | if (status) |
---|
[1121] | 314 | throw utility::GSL_error(std::string("Matrix::mul_elements",status)); |
---|
[716] | 315 | } |
---|
[516] | 316 | |
---|
[716] | 317 | |
---|
[1121] | 318 | void Matrix::resize(size_t r, size_t c, double init_value) |
---|
[808] | 319 | { |
---|
[810] | 320 | delete_allocated_memory(); |
---|
[808] | 321 | |
---|
[1172] | 322 | assert( (r&&c) || (!r&&!c) ); |
---|
| 323 | if (r && c){ |
---|
| 324 | m_ = gsl_matrix_alloc(r,c); |
---|
| 325 | if (!m_) |
---|
| 326 | throw utility::GSL_error("Matrix::Matrix failed to allocate memory"); |
---|
| 327 | all(init_value); |
---|
| 328 | } |
---|
[808] | 329 | |
---|
[810] | 330 | // no need to delete blas_result_ if the number of rows fit, it |
---|
| 331 | // may be useful later. |
---|
| 332 | if (blas_result_ && (blas_result_->size1!=rows())) { |
---|
| 333 | gsl_matrix_free(blas_result_); |
---|
| 334 | blas_result_=NULL; |
---|
| 335 | } |
---|
| 336 | } |
---|
[808] | 337 | |
---|
[810] | 338 | |
---|
[1121] | 339 | size_t Matrix::rows(void) const |
---|
[716] | 340 | { |
---|
[1098] | 341 | return (m_ ? m_->size1 : 0); |
---|
[716] | 342 | } |
---|
| 343 | |
---|
| 344 | |
---|
[1121] | 345 | const VectorConstView Matrix::row_const_view(size_t col) const |
---|
[1009] | 346 | { |
---|
[1027] | 347 | return VectorConstView(*this, col, true); |
---|
[1009] | 348 | } |
---|
| 349 | |
---|
| 350 | |
---|
[1121] | 351 | VectorView Matrix::row_view(size_t row) |
---|
[1009] | 352 | { |
---|
[1015] | 353 | VectorView res(*this, row, true); |
---|
[1009] | 354 | return res; |
---|
| 355 | } |
---|
| 356 | |
---|
| 357 | |
---|
[1121] | 358 | void Matrix::swap_columns(const size_t i, const size_t j) |
---|
[716] | 359 | { |
---|
[792] | 360 | assert(m_); |
---|
[755] | 361 | int status=gsl_matrix_swap_columns(m_, i, j); |
---|
| 362 | if (status) |
---|
[1121] | 363 | throw utility::GSL_error(std::string("Matrix::swap_columns",status)); |
---|
[716] | 364 | } |
---|
| 365 | |
---|
| 366 | |
---|
[1121] | 367 | void Matrix::swap_rowcol(const size_t i, const size_t j) |
---|
[716] | 368 | { |
---|
[792] | 369 | assert(m_); |
---|
[755] | 370 | int status=gsl_matrix_swap_rowcol(m_, i, j); |
---|
| 371 | if (status) |
---|
[1121] | 372 | throw utility::GSL_error(std::string("Matrix::swap_rowcol",status)); |
---|
[716] | 373 | } |
---|
| 374 | |
---|
| 375 | |
---|
[1121] | 376 | void Matrix::swap_rows(const size_t i, const size_t j) |
---|
[716] | 377 | { |
---|
[792] | 378 | assert(m_); |
---|
[755] | 379 | int status=gsl_matrix_swap_rows(m_, i, j); |
---|
| 380 | if (status) |
---|
[1121] | 381 | throw utility::GSL_error(std::string("Matrix::swap_rows",status)); |
---|
[716] | 382 | } |
---|
| 383 | |
---|
| 384 | |
---|
[1121] | 385 | void Matrix::transpose(void) |
---|
[42] | 386 | { |
---|
[792] | 387 | assert(m_); |
---|
[420] | 388 | if (columns()==rows()) |
---|
[754] | 389 | gsl_matrix_transpose(m_); // this never fails |
---|
[420] | 390 | else { |
---|
| 391 | gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows()); |
---|
[753] | 392 | if (!transposed) |
---|
[1121] | 393 | throw utility::GSL_error("Matrix::transpose failed to allocate memory"); |
---|
[754] | 394 | // next line never fails if allocation above succeeded. |
---|
[420] | 395 | gsl_matrix_transpose_memcpy(transposed,m_); |
---|
| 396 | gsl_matrix_free(m_); |
---|
[1098] | 397 | m_ = transposed; |
---|
[762] | 398 | if (blas_result_) { |
---|
| 399 | gsl_matrix_free(blas_result_); |
---|
| 400 | blas_result_=NULL; |
---|
| 401 | } |
---|
[420] | 402 | } |
---|
[42] | 403 | } |
---|
[21] | 404 | |
---|
[12] | 405 | |
---|
[1121] | 406 | double& Matrix::operator()(size_t row, size_t column) |
---|
[716] | 407 | { |
---|
[792] | 408 | assert(m_); |
---|
[813] | 409 | assert(row<rows()); |
---|
| 410 | assert(column<columns()); |
---|
[755] | 411 | double* d=gsl_matrix_ptr(m_, row, column); |
---|
| 412 | if (!d) |
---|
[1121] | 413 | throw utility::GSL_error("Matrix::operator()",GSL_EINVAL); |
---|
[755] | 414 | return *d; |
---|
[716] | 415 | } |
---|
[12] | 416 | |
---|
[716] | 417 | |
---|
[1121] | 418 | const double& Matrix::operator()(size_t row, size_t column) const |
---|
[716] | 419 | { |
---|
[813] | 420 | assert(row<rows()); |
---|
| 421 | assert(column<columns()); |
---|
[1098] | 422 | const double* d=gsl_matrix_const_ptr(m_, row, column); |
---|
[755] | 423 | if (!d) |
---|
[1121] | 424 | throw utility::GSL_error("Matrix::operator()",GSL_EINVAL); |
---|
[755] | 425 | return *d; |
---|
[716] | 426 | } |
---|
| 427 | |
---|
| 428 | |
---|
[1121] | 429 | bool Matrix::operator==(const Matrix& other) const |
---|
[716] | 430 | { |
---|
| 431 | return equal(other); |
---|
| 432 | } |
---|
[421] | 433 | |
---|
[716] | 434 | |
---|
[1121] | 435 | bool Matrix::operator!=(const Matrix& other) const |
---|
[716] | 436 | { |
---|
| 437 | return !equal(other); |
---|
| 438 | } |
---|
| 439 | |
---|
| 440 | |
---|
[1121] | 441 | const Matrix& Matrix::operator=( const Matrix& other ) |
---|
[12] | 442 | { |
---|
[1098] | 443 | if (this!=&other) { |
---|
[1598] | 444 | if (rows()!=other.rows() || columns()!=other.columns()) |
---|
| 445 | resize(other.rows(), other.columns()); |
---|
| 446 | |
---|
| 447 | if (m_) |
---|
| 448 | if (gsl_matrix_memcpy(m_, other.gsl_matrix_p())) { |
---|
| 449 | std::string s="Matrix::create_gsl_matrix_copy dimension mis-match"; |
---|
| 450 | throw utility::GSL_error(s); |
---|
| 451 | } |
---|
[1098] | 452 | } |
---|
[42] | 453 | return *this; |
---|
| 454 | } |
---|
[12] | 455 | |
---|
| 456 | |
---|
[1121] | 457 | const Matrix& Matrix::operator+=(const Matrix& other) |
---|
[703] | 458 | { |
---|
[792] | 459 | assert(m_); |
---|
[1098] | 460 | int status=gsl_matrix_add(m_, other.m_); |
---|
[773] | 461 | if (status) |
---|
[1121] | 462 | throw utility::GSL_error(std::string("Matrix::operator+=", status)); |
---|
[703] | 463 | return *this; |
---|
| 464 | } |
---|
[20] | 465 | |
---|
[703] | 466 | |
---|
[1121] | 467 | const Matrix& Matrix::operator+=(const double d) |
---|
[703] | 468 | { |
---|
[792] | 469 | assert(m_); |
---|
[773] | 470 | gsl_matrix_add_constant(m_, d); |
---|
[703] | 471 | return *this; |
---|
| 472 | } |
---|
| 473 | |
---|
| 474 | |
---|
[1121] | 475 | const Matrix& Matrix::operator-=(const Matrix& other) |
---|
[703] | 476 | { |
---|
[792] | 477 | assert(m_); |
---|
[1098] | 478 | int status=gsl_matrix_sub(m_, other.m_); |
---|
[773] | 479 | if (status) |
---|
[1121] | 480 | throw utility::GSL_error(std::string("Matrix::operator-=", status)); |
---|
[703] | 481 | return *this; |
---|
| 482 | } |
---|
| 483 | |
---|
| 484 | |
---|
[1121] | 485 | const Matrix& Matrix::operator-=(const double d) |
---|
[773] | 486 | { |
---|
[792] | 487 | assert(m_); |
---|
[773] | 488 | gsl_matrix_add_constant(m_, -d); |
---|
| 489 | return *this; |
---|
| 490 | } |
---|
| 491 | |
---|
| 492 | |
---|
[1121] | 493 | const Matrix& Matrix::operator*=(const Matrix& other) |
---|
[42] | 494 | { |
---|
[792] | 495 | assert(m_); |
---|
[796] | 496 | if ( blas_result_ && ((blas_result_->size1!=rows()) || |
---|
| 497 | (blas_result_->size2!=other.columns())) ) { |
---|
| 498 | gsl_matrix_free(blas_result_); |
---|
| 499 | blas_result_=NULL; |
---|
| 500 | } |
---|
[762] | 501 | if (!blas_result_) { |
---|
| 502 | blas_result_ = gsl_matrix_alloc(rows(),other.columns()); |
---|
| 503 | if (!blas_result_) |
---|
[1121] | 504 | throw utility::GSL_error("Matrix::operator*= failed to allocate memory"); |
---|
[762] | 505 | } |
---|
[1098] | 506 | gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, blas_result_); |
---|
[762] | 507 | gsl_matrix* tmp=m_; |
---|
[1098] | 508 | m_ = blas_result_; |
---|
[762] | 509 | blas_result_=tmp; |
---|
[420] | 510 | return *this; |
---|
[42] | 511 | } |
---|
[20] | 512 | |
---|
[12] | 513 | |
---|
[1121] | 514 | const Matrix& Matrix::operator*=(const double d) |
---|
[703] | 515 | { |
---|
[792] | 516 | assert(m_); |
---|
[773] | 517 | gsl_matrix_scale(m_, d); |
---|
[703] | 518 | return *this; |
---|
| 519 | } |
---|
[12] | 520 | |
---|
[703] | 521 | |
---|
[1121] | 522 | bool isnull(const Matrix& other) |
---|
[774] | 523 | { |
---|
[792] | 524 | return gsl_matrix_isnull(other.gsl_matrix_p()); |
---|
[774] | 525 | } |
---|
| 526 | |
---|
| 527 | |
---|
[1121] | 528 | double max(const Matrix& other) |
---|
[774] | 529 | { |
---|
[792] | 530 | return gsl_matrix_max(other.gsl_matrix_p()); |
---|
[774] | 531 | } |
---|
| 532 | |
---|
| 533 | |
---|
[1121] | 534 | double min(const Matrix& other) |
---|
[774] | 535 | { |
---|
[792] | 536 | return gsl_matrix_min(other.gsl_matrix_p()); |
---|
[774] | 537 | } |
---|
| 538 | |
---|
| 539 | |
---|
[1121] | 540 | void minmax_index(const Matrix& other, |
---|
[774] | 541 | std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max) |
---|
| 542 | { |
---|
[792] | 543 | gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second, |
---|
[774] | 544 | &max.first, &max.second); |
---|
| 545 | } |
---|
| 546 | |
---|
| 547 | |
---|
[1121] | 548 | bool nan(const Matrix& templat, Matrix& flag) |
---|
[774] | 549 | { |
---|
| 550 | size_t rows=templat.rows(); |
---|
| 551 | size_t columns=templat.columns(); |
---|
| 552 | if (rows!=flag.rows() && columns!=flag.columns()) |
---|
[1098] | 553 | flag.resize(rows,columns,1.0); |
---|
[774] | 554 | else |
---|
[830] | 555 | flag.all(1.0); |
---|
[774] | 556 | bool nan=false; |
---|
| 557 | for (size_t i=0; i<rows; i++) |
---|
| 558 | for (size_t j=0; j<columns; j++) |
---|
| 559 | if (std::isnan(templat(i,j))) { |
---|
| 560 | flag(i,j)=0; |
---|
| 561 | nan=true; |
---|
| 562 | } |
---|
| 563 | return nan; |
---|
| 564 | } |
---|
| 565 | |
---|
| 566 | |
---|
[1121] | 567 | void swap(Matrix& a, Matrix& b) |
---|
[774] | 568 | { |
---|
[794] | 569 | assert(a.gsl_matrix_p()); assert(b.gsl_matrix_p()); |
---|
[774] | 570 | int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p()); |
---|
| 571 | if (status) |
---|
[1121] | 572 | throw utility::GSL_error(std::string("swap(Matrix&,Matrix&)",status)); |
---|
[774] | 573 | } |
---|
| 574 | |
---|
| 575 | |
---|
[1121] | 576 | std::ostream& operator<<(std::ostream& s, const Matrix& m) |
---|
[42] | 577 | { |
---|
[420] | 578 | s.setf(std::ios::dec); |
---|
[42] | 579 | s.precision(12); |
---|
[420] | 580 | for(size_t i=0, j=0; i<m.rows(); i++) |
---|
[273] | 581 | for (j=0; j<m.columns(); j++) { |
---|
| 582 | s << m(i,j); |
---|
| 583 | if (j<m.columns()-1) |
---|
[593] | 584 | s << s.fill(); |
---|
[273] | 585 | else if (i<m.rows()-1) |
---|
| 586 | s << "\n"; |
---|
| 587 | } |
---|
[420] | 588 | return s; |
---|
[42] | 589 | } |
---|
[12] | 590 | |
---|
[737] | 591 | |
---|
[1121] | 592 | Vector operator*(const Matrix& m, const VectorBase& v) |
---|
[737] | 593 | { |
---|
[1120] | 594 | utility::Vector res(m.rows()); |
---|
[737] | 595 | for (size_t i=0; i<res.size(); ++i) |
---|
[1027] | 596 | res(i) = VectorConstView(m,i) * v; |
---|
[737] | 597 | return res; |
---|
| 598 | } |
---|
| 599 | |
---|
| 600 | |
---|
[1121] | 601 | Vector operator*(const VectorBase& v, const Matrix& m) |
---|
[737] | 602 | { |
---|
[1120] | 603 | utility::Vector res(m.columns()); |
---|
[737] | 604 | for (size_t i=0; i<res.size(); ++i) |
---|
[1027] | 605 | res(i) = v * VectorConstView(m,i,false); |
---|
[737] | 606 | return res; |
---|
| 607 | } |
---|
| 608 | |
---|
[680] | 609 | }}} // of namespace utility, yat and thep |
---|