source: trunk/yat/utility/Matrix.cc @ 1598

Last change on this file since 1598 was 1598, checked in by Peter, 13 years ago

Fixed some issues related to empty (zero sized) Matrix.

  • Matrix(size_t, size_t) now accepts 0 as argument
  • Matrix(const& Matrix) now works when passed Matrix is empty
  • assignment now works when rhs is empty

related to ticket:343

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.2 KB
RevLine 
[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
39namespace theplu {
[680]40namespace yat {
[616]41namespace 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
Note: See TracBrowser for help on using the repository browser.