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

Last change on this file since 1172 was 1172, checked in by Peter, 14 years ago

allowing resize to zero size. In Matrix if resize is called with only one of columns or rows set to zero a zero by zero matrix is created. I think we can live with that. Though an assertion is called because it is strange.

  • 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 1172 2008-02-27 14:32:15Z 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, Markus Ringnér, Peter Johansson
7  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://trac.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    u_int nof_columns=0;
79    u_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(utility::is_double(element)) {
101          v.push_back(atof(element.c_str()));
102        }
103        else if (!element.size() || utility::is_nan(element)) {
104          v.push_back(std::numeric_limits<double>::quiet_NaN());
105        }
106        else {
107          std::stringstream ss("Warning: '");
108          ss << element << "' is not accepted as a matrix element.";
109          throw IO_error(ss.str());
110        }
111      }           
112      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
113          v.push_back(std::numeric_limits<double>::quiet_NaN());
114      if (!nof_columns)
115        nof_columns=v.size();
116      else if (v.size()!=nof_columns) {
117        std::ostringstream s;
118        s << "Matrix::Matrix(std::istream&, char) data file error: "
119          << "line " << nof_rows << " has " << v.size()
120          << " columns; expected " << nof_columns << " columns.";
121        throw utility::IO_error(s.str());
122      }
123      data_matrix.push_back(v);
124    }
125
126    // manipulate the state of the stream to be good
127    is.clear(std::ios::goodbit);
128   
129    // if stream was empty, create nothing
130    if (!nof_columns || !nof_rows)
131      return;
132
133    // convert the data to a gsl matrix
134    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
135    if (!m_)
136      throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
137
138    // if gsl error handler disabled, out of bounds index will not
139    // abort the program.
140    for(u_int i=0;i<nof_rows;i++)
141      for(u_int j=0;j<nof_columns;j++)
142        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
143  }
144
145
146  Matrix::~Matrix(void)
147  {
148    delete_allocated_memory();
149    if (blas_result_)
150      gsl_matrix_free(blas_result_);
151    blas_result_=NULL;
152  }
153
154
155  void Matrix::all(const double value)
156  {
157    assert(m_);
158    gsl_matrix_set_all(m_, value);
159  }
160
161
162  Matrix::iterator Matrix::begin(void)
163  {
164    return iterator(&(*this)(0,0), 1);
165  }
166
167
168  Matrix::const_iterator Matrix::begin(void) const
169  {
170    return const_iterator(&(*this)(0,0), 1);
171  }
172
173
174  Matrix::column_iterator Matrix::begin_column(size_t i)
175  {
176    return iterator(&(*this)(0,i), this->columns());
177  }
178
179
180  Matrix::const_column_iterator Matrix::begin_column(size_t i) const
181  {
182    return const_iterator(&(*this)(0,i), this->columns());
183  }
184
185
186  Matrix::row_iterator Matrix::begin_row(size_t i)
187  {
188    return iterator(&(*this)(i,0), 1);
189  }
190
191
192  Matrix::const_row_iterator Matrix::begin_row(size_t i) const
193  {
194    return const_iterator(&(*this)(i,0), 1);
195  }
196
197
198  VectorView Matrix::column_view(size_t col)
199  {
200    VectorView res(*this, col, false);
201    return res;
202  }
203
204
205  const VectorConstView Matrix::column_const_view(size_t col) const
206  {
207    return VectorConstView(*this, col, false);
208  }
209
210
211  size_t Matrix::columns(void) const
212  {
213    return (m_ ? m_->size2 : 0);
214  }
215
216
217  gsl_matrix* Matrix::create_gsl_matrix_copy(void) const
218  {
219    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
220    if (!m)
221      throw utility::GSL_error("Matrix::create_gsl_matrix_copy failed to allocate memory");
222    if (gsl_matrix_memcpy(m,m_))
223      throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
224    return m;
225  }
226
227
228  void Matrix::delete_allocated_memory(void)
229  {
230    if (m_)
231      gsl_matrix_free(m_);
232    m_=NULL;
233  }
234
235
236  void Matrix::div(const Matrix& other)
237  {
238    assert(m_);
239    int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
240    if (status)
241      throw utility::GSL_error(std::string("matrix::div_elements",status));
242  }
243
244
245  Matrix::iterator Matrix::end(void)
246  {
247    return iterator(&(*this)(0,0)+rows()*columns(), 1);
248  }
249
250
251  Matrix::const_iterator Matrix::end(void) const
252  {
253    return const_iterator(&(*this)(0,0)+rows()*columns(), 1);
254  }
255
256
257  Matrix::column_iterator Matrix::end_column(size_t i)
258  {
259    return column_iterator(&(*this)(0,i)+rows()*columns(), this->columns());
260  }
261
262
263  Matrix::const_column_iterator Matrix::end_column(size_t i) const
264  {
265    return const_column_iterator(&(*this)(0,i)+rows()*columns(),this->columns());
266  }
267
268
269  Matrix::row_iterator Matrix::end_row(size_t i)
270  {
271    return row_iterator(&(*this)(i,0)+columns(), 1);
272  }
273
274
275  Matrix::const_row_iterator Matrix::end_row(size_t i) const
276  {
277    return const_row_iterator(&(*this)(i,0)+columns(), 1);
278  }
279
280
281  bool Matrix::equal(const Matrix& other, const double d) const
282  {
283    if (this==&other)
284      return true;
285    if (columns()!=other.columns() || rows()!=other.rows())
286      return false;
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 ||
291            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
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    assert(other.m_);
444    if (this!=&other) {
445      if ( !m_ || (other.m_->size1!=m_->size1) || (other.m_->size2!=m_->size2) )
446        resize(other.m_->size1,other.m_->size2);
447      if (gsl_matrix_memcpy(m_, other.gsl_matrix_p()))
448        throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
449    }
450    return *this;
451  }
452
453
454  const Matrix& Matrix::operator+=(const Matrix& other)
455  {
456    assert(m_);
457    int status=gsl_matrix_add(m_, other.m_);
458    if (status)
459      throw utility::GSL_error(std::string("Matrix::operator+=", status));
460    return *this;
461  }
462
463
464  const Matrix& Matrix::operator+=(const double d)
465  {
466    assert(m_);
467    gsl_matrix_add_constant(m_, d);
468    return *this;
469  }
470
471
472  const Matrix& Matrix::operator-=(const Matrix& other)
473  {
474    assert(m_);
475    int status=gsl_matrix_sub(m_, other.m_);
476    if (status)
477      throw utility::GSL_error(std::string("Matrix::operator-=", status));
478    return *this;
479  }
480
481
482  const Matrix& Matrix::operator-=(const double d)
483  {
484    assert(m_);
485    gsl_matrix_add_constant(m_, -d);
486    return *this;
487  }
488
489
490  const Matrix& Matrix::operator*=(const Matrix& other)
491  {
492    assert(m_);
493    if ( blas_result_ && ((blas_result_->size1!=rows()) ||
494                          (blas_result_->size2!=other.columns())) ) {
495      gsl_matrix_free(blas_result_);
496      blas_result_=NULL;
497    }
498    if (!blas_result_) {
499      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
500      if (!blas_result_)
501        throw utility::GSL_error("Matrix::operator*= failed to allocate memory");
502    }
503    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, blas_result_);
504    gsl_matrix* tmp=m_;
505    m_ = blas_result_;
506    blas_result_=tmp;
507    return *this;
508  }
509
510
511  const Matrix& Matrix::operator*=(const double d)
512  {
513    assert(m_);
514    gsl_matrix_scale(m_, d);
515    return *this;
516  }
517
518
519  bool isnull(const Matrix& other)
520  {
521    return gsl_matrix_isnull(other.gsl_matrix_p());
522  }
523
524
525  double max(const Matrix& other)
526  {
527    return gsl_matrix_max(other.gsl_matrix_p());
528  }
529
530
531  double min(const Matrix& other)
532  {
533    return gsl_matrix_min(other.gsl_matrix_p());
534  }
535
536
537  void minmax_index(const Matrix& other,
538                    std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
539  {
540    gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,
541                            &max.first, &max.second);
542  }
543
544
545  bool nan(const Matrix& templat, Matrix& flag)
546  {
547    size_t rows=templat.rows();
548    size_t columns=templat.columns();
549    if (rows!=flag.rows() && columns!=flag.columns())
550      flag.resize(rows,columns,1.0);
551    else
552      flag.all(1.0);
553    bool nan=false;
554    for (size_t i=0; i<rows; i++)
555      for (size_t j=0; j<columns; j++) 
556        if (std::isnan(templat(i,j))) {
557          flag(i,j)=0;
558          nan=true;
559        }
560    return nan;
561  }
562
563
564  void swap(Matrix& a, Matrix& b)
565  {
566    assert(a.gsl_matrix_p()); assert(b.gsl_matrix_p());
567    int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
568    if (status)
569      throw utility::GSL_error(std::string("swap(Matrix&,Matrix&)",status));
570  }
571
572
573  std::ostream& operator<<(std::ostream& s, const Matrix& m)
574  {
575    s.setf(std::ios::dec);
576    s.precision(12);
577    for(size_t i=0, j=0; i<m.rows(); i++)
578      for (j=0; j<m.columns(); j++) {
579        s << m(i,j);
580        if (j<m.columns()-1)
581          s << s.fill();
582        else if (i<m.rows()-1)
583          s << "\n";
584      }
585    return s;
586  }
587
588
589  Vector operator*(const Matrix& m, const VectorBase& v)
590  {
591    utility::Vector res(m.rows());
592    for (size_t i=0; i<res.size(); ++i)
593      res(i) = VectorConstView(m,i) * v;
594    return res;
595  }
596
597
598  Vector operator*(const VectorBase& v, const Matrix& m)
599  {
600    utility::Vector res(m.columns());
601    for (size_t i=0; i<res.size(); ++i)
602      res(i) = v * VectorConstView(m,i,false);
603    return res;
604  }
605
606}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.