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

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

merge patch release 0.4.2 to trunk. Delta 0.4.2-0.4.1

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.3 KB
Line 
1// $Id: Matrix.cc 1437 2008-08-25 17:55:00Z 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 2 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 this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
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"
33
34#include <cassert>
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_(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 istream
68  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 (dynamically
76    // 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 lines
83      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        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      data_matrix.push_back(v);
126    }
127
128    // manipulate the state of the stream to be good
129    is.clear(std::ios::goodbit);
130   
131    // if stream was empty, create nothing
132    if (!nof_columns || !nof_rows)
133      return;
134
135    // convert the data to a gsl matrix
136    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 not
141    // 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) const
171  {
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) const
183  {
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) const
195  {
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) const
208  {
209    return VectorConstView(*this, col, false);
210  }
211
212
213  size_t Matrix::columns(void) const
214  {
215    return (m_ ? m_->size2 : 0);
216  }
217
218
219  gsl_matrix* Matrix::create_gsl_matrix_copy(void) const
220  {
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 mis-match");
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) const
254  {
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) const
266  {
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) const
278  {
279    return const_row_iterator(&(*this)(i,0)+columns(), 1);
280  }
281
282
283  bool Matrix::equal(const Matrix& other, const double d) const
284  {
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 detection
292        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) const
300  {
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, it
333    // 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) const
342  {
343    return (m_ ? m_->size1 : 0);
344  }
345
346
347  const VectorConstView Matrix::row_const_view(size_t col) const
348  {
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 fails
392    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) const
421  {
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) const
432  {
433    return equal(other);
434  }
435
436
437  bool Matrix::operator!=(const Matrix& other) const
438  {
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 mis-match");
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    else
554      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}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.