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

Last change on this file since 1786 was 1786, checked in by Peter, 12 years ago

fixes #421

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.0 KB
Line 
1// $Id: Matrix.cc 1786 2009-02-09 20:35:25Z peter $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
5  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
7  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://dev.thep.lu.se/yat
10
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
13  published by the Free Software Foundation; either version 3 of the
14  License, or (at your option) any later version.
15
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
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with yat. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include "Matrix.h"
26#include "stl_utility.h"
27#include "Vector.h"
28#include "VectorBase.h"
29#include "VectorConstView.h"
30#include "VectorView.h"
31#include "utility.h"
32
33#include <cassert>
34#include <climits>
35#include <cmath>
36#include <sstream>
37#include <vector>
38
39#include <gsl/gsl_blas.h>
40
41namespace theplu {
42namespace yat {
43namespace utility {
44
45
46  Matrix::Matrix(void)
47    : blas_result_(NULL), m_(NULL)
48  {
49  }
50
51
52  Matrix::Matrix(const size_t& r, const size_t& c, double init_value)
53    : blas_result_(NULL), m_(NULL)
54  {
55    resize(r,c,init_value);
56  }
57
58
59  Matrix::Matrix(const Matrix& o)
60    : blas_result_(NULL), m_(o.create_gsl_matrix_copy())
61  {
62  }
63
64
65  // Constructor that gets data from istream
66  Matrix::Matrix(std::istream& is, char sep) 
67    throw (utility::IO_error,std::exception)
68    : blas_result_(NULL)
69  {
70    if (!is.good())
71      throw utility::IO_error("Matrix: istream is not good");
72
73    // read the data file and store in stl vectors (dynamically
74    // expandable)
75    std::vector<std::vector<double> > data_matrix;
76    unsigned int nof_columns=0;
77    unsigned int nof_rows = 0;
78    std::string line;
79    while(getline(is, line, '\n')){
80      // Ignoring empty lines
81      if (!line.size()) {
82        continue;
83      }
84      nof_rows++;
85      data_matrix.resize(data_matrix.size()+1);
86      std::vector<double>& v=data_matrix.back();
87      v.reserve(nof_columns);
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        else
96          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: " << element
109               << " 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 separator
115          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    }
126
127    // manipulate the state of the stream to be good
128    is.clear(std::ios::goodbit);
129   
130    // if stream was empty, create nothing
131    if (!nof_columns || !nof_rows)
132      return;
133
134    // convert the data to a gsl matrix
135    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
136    if (!m_)
137      throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
138
139    // if gsl error handler disabled, out of bounds index will not
140    // abort the program.
141    for(size_t i=0;i<nof_rows;i++)
142      for(size_t j=0;j<nof_columns;j++)
143        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
144  }
145
146
147  Matrix::~Matrix(void)
148  {
149    delete_allocated_memory();
150    if (blas_result_)
151      gsl_matrix_free(blas_result_);
152    blas_result_=NULL;
153  }
154
155
156  void Matrix::all(const double value)
157  {
158    assert(m_);
159    gsl_matrix_set_all(m_, value);
160  }
161
162
163  Matrix::iterator Matrix::begin(void)
164  {
165    return iterator(&(*this)(0,0), 1);
166  }
167
168
169  Matrix::const_iterator Matrix::begin(void) const
170  {
171    return const_iterator(&(*this)(0,0), 1);
172  }
173
174
175  Matrix::column_iterator Matrix::begin_column(size_t i)
176  {
177    return column_iterator(&(*this)(0,i), this->columns());
178  }
179
180
181  Matrix::const_column_iterator Matrix::begin_column(size_t i) const
182  {
183    return const_column_iterator(&(*this)(0,i), this->columns());
184  }
185
186
187  Matrix::row_iterator Matrix::begin_row(size_t i)
188  {
189    return row_iterator(&(*this)(i,0), 1);
190  }
191
192
193  Matrix::const_row_iterator Matrix::begin_row(size_t i) const
194  {
195    return const_row_iterator(&(*this)(i,0), 1);
196  }
197
198
199  VectorView Matrix::column_view(size_t col)
200  {
201    VectorView res(*this, col, false);
202    return res;
203  }
204
205
206  const VectorConstView Matrix::column_const_view(size_t col) const
207  {
208    return VectorConstView(*this, col, false);
209  }
210
211
212  size_t Matrix::columns(void) const
213  {
214    return (m_ ? m_->size2 : 0);
215  }
216
217
218  gsl_matrix* Matrix::create_gsl_matrix_copy(void) const
219  {
220    if (!m_)
221      return NULL;
222    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
223    if (!m)
224      throw utility::GSL_error("Matrix::create_gsl_matrix_copy failed to allocate memory");
225    if (gsl_matrix_memcpy(m,m_))
226      throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
227    return m;
228  }
229
230
231  void Matrix::delete_allocated_memory(void)
232  {
233    if (m_)
234      gsl_matrix_free(m_);
235    m_=NULL;
236  }
237
238
239  void Matrix::div(const Matrix& other)
240  {
241    assert(m_);
242    int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
243    if (status)
244      throw utility::GSL_error(std::string("matrix::div_elements",status));
245  }
246
247
248  Matrix::iterator Matrix::end(void)
249  {
250    return iterator(&(*this)(0,0)+rows()*columns(), 1);
251  }
252
253
254  Matrix::const_iterator Matrix::end(void) const
255  {
256    return const_iterator(&(*this)(0,0)+rows()*columns(), 1);
257  }
258
259
260  Matrix::column_iterator Matrix::end_column(size_t i)
261  {
262    return column_iterator(&(*this)(0,i)+rows()*columns(), this->columns());
263  }
264
265
266  Matrix::const_column_iterator Matrix::end_column(size_t i) const
267  {
268    return const_column_iterator(&(*this)(0,i)+rows()*columns(),this->columns());
269  }
270
271
272  Matrix::row_iterator Matrix::end_row(size_t i)
273  {
274    return row_iterator(&(*this)(i,0)+columns(), 1);
275  }
276
277
278  Matrix::const_row_iterator Matrix::end_row(size_t i) const
279  {
280    return const_row_iterator(&(*this)(i,0)+columns(), 1);
281  }
282
283
284  bool Matrix::equal(const Matrix& other, const double d) const
285  {
286    if (this==&other)
287      return true;
288    if (columns()!=other.columns() || rows()!=other.rows())
289      return false;
290    for (size_t i=0; i<rows(); i++)
291      if (!row_const_view(i).equal(other.row_const_view(i), d))
292        return false;
293    return true;
294  }
295
296
297  const gsl_matrix* Matrix::gsl_matrix_p(void) const
298  {
299    return m_;
300  }
301
302
303  gsl_matrix* Matrix::gsl_matrix_p(void)
304  {
305    return m_;
306  }
307
308
309  void Matrix::mul(const Matrix& other)
310  {
311    assert(m_);
312    int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());
313    if (status)
314      throw utility::GSL_error(std::string("Matrix::mul_elements",status));
315  }
316
317
318  void Matrix::resize(size_t r, size_t c, double init_value)
319  {
320    delete_allocated_memory();
321
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    }
329
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  }
337
338
339  size_t Matrix::rows(void) const
340  {
341    return (m_ ? m_->size1 : 0);
342  }
343
344
345  const VectorConstView Matrix::row_const_view(size_t col) const
346  {
347    return VectorConstView(*this, col, true);
348  }
349
350
351  VectorView Matrix::row_view(size_t row)
352  {
353    VectorView res(*this, row, true);
354    return res;
355  }
356
357
358  void Matrix::swap_columns(const size_t i, const size_t j)
359  {
360    assert(m_);
361    int status=gsl_matrix_swap_columns(m_, i, j);
362    if (status)
363      throw utility::GSL_error(std::string("Matrix::swap_columns",status));
364  }
365
366
367  void Matrix::swap_rowcol(const size_t i, const size_t j)
368  {
369    assert(m_);
370    int status=gsl_matrix_swap_rowcol(m_, i, j);
371    if (status)
372      throw utility::GSL_error(std::string("Matrix::swap_rowcol",status));
373  }
374
375
376  void Matrix::swap_rows(const size_t i, const size_t j)
377  {
378    assert(m_);
379    int status=gsl_matrix_swap_rows(m_, i, j);
380    if (status)
381      throw utility::GSL_error(std::string("Matrix::swap_rows",status));
382  }
383
384
385  void Matrix::transpose(void)
386  {
387    assert(m_);
388    if (columns()==rows())
389      gsl_matrix_transpose(m_); // this never fails
390    else {
391      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
392      if (!transposed)
393        throw utility::GSL_error("Matrix::transpose failed to allocate memory");
394      // next line never fails if allocation above succeeded.
395      gsl_matrix_transpose_memcpy(transposed,m_);
396      gsl_matrix_free(m_);
397      m_ = transposed;
398      if (blas_result_) {
399        gsl_matrix_free(blas_result_);
400        blas_result_=NULL;
401      }
402    }
403  }
404
405
406  double& Matrix::operator()(size_t row, size_t column)
407  {
408    assert(m_);
409    assert(row<rows());
410    assert(column<columns());
411    double* d=gsl_matrix_ptr(m_, row, column);
412    if (!d)
413      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
414    return *d;
415  }
416
417
418  const double& Matrix::operator()(size_t row, size_t column) const
419  {
420    assert(row<rows());
421    assert(column<columns());
422    const double* d=gsl_matrix_const_ptr(m_, row, column);
423    if (!d)
424      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
425    return *d;
426  }
427
428
429  bool Matrix::operator==(const Matrix& other) const
430  {
431    return equal(other);
432  }
433
434
435  bool Matrix::operator!=(const Matrix& other) const
436  {
437    return !equal(other);
438  }
439
440
441  const Matrix& Matrix::operator=( const Matrix& other )
442  {
443    if (this!=&other) {
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::operator=  assignment failed";
450          throw utility::GSL_error(s);
451        }
452    }
453    return *this;
454  }
455
456
457  const Matrix& Matrix::operator+=(const Matrix& other)
458  {
459    assert(m_);
460    int status=gsl_matrix_add(m_, other.m_);
461    if (status)
462      throw utility::GSL_error(std::string("Matrix::operator+=", status));
463    return *this;
464  }
465
466
467  const Matrix& Matrix::operator+=(const double d)
468  {
469    assert(m_);
470    gsl_matrix_add_constant(m_, d);
471    return *this;
472  }
473
474
475  const Matrix& Matrix::operator-=(const Matrix& other)
476  {
477    assert(m_);
478    int status=gsl_matrix_sub(m_, other.m_);
479    if (status)
480      throw utility::GSL_error(std::string("Matrix::operator-=", status));
481    return *this;
482  }
483
484
485  const Matrix& Matrix::operator-=(const double d)
486  {
487    assert(m_);
488    gsl_matrix_add_constant(m_, -d);
489    return *this;
490  }
491
492
493  const Matrix& Matrix::operator*=(const Matrix& other)
494  {
495    assert(m_);
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    }
501    if (!blas_result_) {
502      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
503      if (!blas_result_)
504        throw utility::GSL_error("Matrix::operator*= failed to allocate memory");
505    }
506    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, blas_result_);
507    gsl_matrix* tmp=m_;
508    m_ = blas_result_;
509    blas_result_=tmp;
510    return *this;
511  }
512
513
514  const Matrix& Matrix::operator*=(const double d)
515  {
516    assert(m_);
517    gsl_matrix_scale(m_, d);
518    return *this;
519  }
520
521
522  bool isnull(const Matrix& other)
523  {
524    return gsl_matrix_isnull(other.gsl_matrix_p());
525  }
526
527
528  double max(const Matrix& other)
529  {
530    return gsl_matrix_max(other.gsl_matrix_p());
531  }
532
533
534  double min(const Matrix& other)
535  {
536    return gsl_matrix_min(other.gsl_matrix_p());
537  }
538
539
540  void minmax_index(const Matrix& other,
541                    std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
542  {
543    gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,
544                            &max.first, &max.second);
545  }
546
547
548  bool nan(const Matrix& templat, Matrix& flag)
549  {
550    if (templat.rows()!=flag.rows() || templat.columns()!=flag.columns())
551      flag.resize(templat.rows(),templat.columns(),1.0);
552    return BinaryWeight()(templat.begin(), templat.end(), flag.begin());
553  }
554
555
556  void swap(Matrix& a, Matrix& b)
557  {
558    assert(a.gsl_matrix_p()); assert(b.gsl_matrix_p());
559    int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
560    if (status)
561      throw utility::GSL_error(std::string("swap(Matrix&,Matrix&)",status));
562  }
563
564
565  std::ostream& operator<<(std::ostream& s, const Matrix& m)
566  {
567    s.setf(std::ios::dec);
568    s.precision(12);
569    for(size_t i=0, j=0; i<m.rows(); i++)
570      for (j=0; j<m.columns(); j++) {
571        s << m(i,j);
572        if (j<m.columns()-1)
573          s << s.fill();
574        else if (i<m.rows()-1)
575          s << "\n";
576      }
577    return s;
578  }
579
580
581  Vector operator*(const Matrix& m, const VectorBase& v)
582  {
583    utility::Vector res(m.rows());
584    for (size_t i=0; i<res.size(); ++i)
585      res(i) = VectorConstView(m,i) * v;
586    return res;
587  }
588
589
590  Vector operator*(const VectorBase& v, const Matrix& m)
591  {
592    utility::Vector res(m.columns());
593    for (size_t i=0; i<res.size(); ++i)
594      res(i) = v * VectorConstView(m,i,false);
595    return res;
596  }
597
598}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.