source: trunk/yat/utility/matrix.cc @ 788

Last change on this file since 788 was 788, checked in by Jari Häkkinen, 16 years ago

Improved comparison operator and equal function doxygen doc. Added for self comparison in equal.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.1 KB
Line 
1// $Id: matrix.cc 788 2007-03-10 19:59:51Z jari $
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 Jari Häkkinen
8
9  This file is part of the yat library, http://lev.thep.lu.se/trac/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 "utility.h"
30
31#include <cmath>
32#include <sstream>
33#include <vector>
34
35#include <gsl/gsl_blas.h>
36
37namespace theplu {
38namespace yat {
39namespace utility {
40
41
42  matrix::matrix(void)
43    : blas_result_(NULL), m_(NULL), view_(NULL)
44  {
45  }
46
47
48  matrix::matrix(const size_t& r, const size_t& c, double init_value)
49    : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL)
50  {
51    if (!m_)
52      throw utility::GSL_error("matrix::matrix failed to allocate memory");
53    set_all(init_value);
54  }
55
56
57  matrix::matrix(const matrix& o)
58    : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL)
59  {
60  }
61
62
63  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
64                 size_t n_row, size_t n_column)
65    : blas_result_(NULL)
66  {
67    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
68                                                     offset_row,offset_column,
69                                                     n_row,n_column));
70    if (!view_)
71      throw utility::GSL_error("matrix::matrix failed to setup view");
72    m_ = &(view_->matrix);
73  }
74
75
76  // Constructor that gets data from istream
77  matrix::matrix(std::istream& is, char sep) 
78    throw (utility::IO_error,std::exception)
79    : blas_result_(NULL), view_(NULL)
80  {
81    // read the data file and store in stl vectors (dynamically
82    // expandable)
83    std::vector<std::vector<double> > data_matrix;
84    u_int nof_columns=0;
85    u_int nof_rows = 0;
86    std::string line;
87    while(getline(is, line, '\n')){
88      // Ignoring empty lines
89      if (!line.size()) {
90        continue;
91      }
92      nof_rows++;
93      std::vector<double> v;
94      std::string element;
95      std::stringstream ss(line);
96     
97      bool ok=true;
98      while(ok) {
99        if(sep=='\0')
100          ok=(ss>>element);
101        else
102          ok=getline(ss, element, sep);
103        if(!ok)
104          break;
105       
106        if(utility::is_double(element)) {
107          v.push_back(atof(element.c_str()));
108        }
109        else if (!element.size() || utility::is_nan(element)) {
110          v.push_back(std::numeric_limits<double>::quiet_NaN());
111        }
112        else {
113          std::stringstream ss("Warning: '");
114          ss << element << "' is not accepted as a matrix element.";
115          throw IO_error(ss.str());
116        }
117      }           
118      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
119          v.push_back(std::numeric_limits<double>::quiet_NaN());
120      if (!nof_columns)
121        nof_columns=v.size();
122      else if (v.size()!=nof_columns) {
123        std::ostringstream s;
124        s << "matrix::matrix(std::istream&, char) data file error: "
125          << "line " << nof_rows << " has " << v.size()
126          << " columns; expected " << nof_columns << " columns.";
127        throw utility::IO_error(s.str());
128      }
129      data_matrix.push_back(v);
130    }
131
132    // manipulate the state of the stream to be good
133    is.clear(std::ios::goodbit);
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(u_int i=0;i<nof_rows;i++)
142      for(u_int 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    if (view_)
150      delete view_;
151    else {
152      if (blas_result_)
153        gsl_matrix_free(blas_result_);
154      if (m_)
155        gsl_matrix_free(m_);
156    }
157    blas_result_=NULL;
158    m_=NULL;
159  }
160
161
162  size_t matrix::columns(void) const
163  {
164    if (!m_)
165      return 0;
166    return m_->size2;
167  }
168
169
170  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
171  {
172    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
173    if (!m)
174      throw utility::GSL_error("matrix::create_gsl_matrix_copy failed to allocate memory");
175    if (gsl_matrix_memcpy(m,m_))
176      throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mis-match");
177    return m;
178  }
179
180
181  void matrix::div_elements(const matrix& b)
182  {
183    int status=gsl_matrix_div_elements(m_, b.m_);
184    if (status)
185      throw utility::GSL_error(std::string("matrix::div_elements",status));
186  }
187
188
189  bool matrix::equal(const matrix& other, const double d) const
190  {
191    if (this==&other)
192      return true;
193    if (columns()!=other.columns() || rows()!=other.rows())
194      return false;
195    for (size_t i=0; i<rows(); i++)
196      for (size_t j=0; j<columns(); j++)
197        // The two last condition checks are needed for NaN detection
198        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
199            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
200          return false;
201    return true;
202  }
203
204
205  const gsl_matrix* matrix::gsl_matrix_p(void) const
206  {
207    return m_;
208  }
209
210
211  gsl_matrix* matrix::gsl_matrix_p(void)
212  {
213    return m_;
214  }
215
216
217  bool matrix::isview(void) const
218  {
219    return view_;
220  }
221
222
223  void matrix::mul_elements(const matrix& b)
224  {
225    int status=gsl_matrix_mul_elements(m_, b.m_);
226    if (status)
227      throw utility::GSL_error(std::string("matrix::mul_elements",status));
228   
229  }
230
231
232  size_t matrix::rows(void) const
233  {
234    if (!m_)
235      return 0;
236    return m_->size1;
237  }
238
239
240  void matrix::set(const matrix& mat)
241  {
242    if (gsl_matrix_memcpy(m_, mat.m_))
243      throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mis-match");
244  }
245
246
247  void matrix::set_all(const double value)
248  {
249    gsl_matrix_set_all(m_, value);
250  }
251
252
253  void matrix::set_column(const size_t column, const vector& vec)
254  {
255    int status=gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
256    if (status)
257      throw utility::GSL_error(std::string("matrix::set_column",status));
258  }
259
260
261  void matrix::set_row(const size_t row, const vector& vec)
262  {
263    int status=gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
264    if (status)
265      throw utility::GSL_error(std::string("matrix::set_row",status));
266  }
267
268
269  void matrix::swap_columns(const size_t i, const size_t j)
270  {
271    int status=gsl_matrix_swap_columns(m_, i, j);
272    if (status)
273      throw utility::GSL_error(std::string("matrix::swap_columns",status));
274  }
275
276
277  void matrix::swap_rowcol(const size_t i, const size_t j)
278  {
279    int status=gsl_matrix_swap_rowcol(m_, i, j);
280    if (status)
281      throw utility::GSL_error(std::string("matrix::swap_rowcol",status));
282  }
283
284
285  void matrix::swap_rows(const size_t i, const size_t j)
286  {
287    int status=gsl_matrix_swap_rows(m_, i, j);
288    if (status)
289      throw utility::GSL_error(std::string("matrix::swap_rows",status));
290  }
291
292
293  void matrix::transpose(void)
294  {
295    if (columns()==rows())
296      gsl_matrix_transpose(m_); // this never fails
297    else {
298      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
299      if (!transposed)
300        throw utility::GSL_error("matrix::transpose failed to allocate memory");
301      // next line never fails if allocation above succeeded.
302      gsl_matrix_transpose_memcpy(transposed,m_);
303      gsl_matrix_free(m_);
304      m_=transposed;
305      if (blas_result_) {
306        gsl_matrix_free(blas_result_);
307        blas_result_=NULL;
308      }
309    }
310  }
311
312
313  double& matrix::operator()(size_t row, size_t column)
314  {
315    double* d=gsl_matrix_ptr(m_, row, column);
316    if (!d)
317      throw utility::GSL_error("matrix::operator()",GSL_EINVAL);
318    return *d;
319  }
320
321
322  const double& matrix::operator()(size_t row, size_t column) const
323  {
324    const double* d=gsl_matrix_const_ptr(m_, row, column);
325    if (!d)
326      throw utility::GSL_error("matrix::operator()",GSL_EINVAL);
327    return *d;
328  }
329
330
331  bool matrix::operator==(const matrix& other) const
332  {
333    return equal(other);
334  }
335
336
337  bool matrix::operator!=(const matrix& other) const
338  {
339    return !equal(other);
340  }
341
342
343  const matrix& matrix::operator=( const matrix& other )
344  {
345    if ( this != &other ) {
346      if (view_)
347        delete view_;
348      else if (m_)
349        gsl_matrix_free(m_);
350      m_ = other.create_gsl_matrix_copy();
351    }
352    return *this;
353  }
354
355
356  const matrix& matrix::operator+=(const matrix& m)
357  {
358    int status=gsl_matrix_add(m_, m.m_);
359    if (status)
360      throw utility::GSL_error(std::string("matrix::operator+=", status));
361    return *this;
362  }
363
364
365  const matrix& matrix::operator+=(const double d)
366  {
367    gsl_matrix_add_constant(m_, d);
368    return *this;
369  }
370
371
372  const matrix& matrix::operator-=(const matrix& m)
373  {
374    int status=gsl_matrix_sub(m_, m.m_);
375    if (status)
376      throw utility::GSL_error(std::string("matrix::operator-=", status));
377    return *this;
378  }
379
380
381  const matrix& matrix::operator-=(const double d)
382  {
383    gsl_matrix_add_constant(m_, -d);
384    return *this;
385  }
386
387
388  const matrix& matrix::operator*=(const matrix& other)
389  {
390    if (!blas_result_) {
391      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
392      if (!blas_result_)
393        throw utility::GSL_error("matrix::operator*= failed to allocate memory");
394    }
395    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0,
396                   blas_result_);
397    gsl_matrix* tmp=m_;
398    m_=blas_result_;
399    blas_result_=tmp;
400    return *this;
401  }
402
403
404  const matrix& matrix::operator*=(const double d)
405  {
406    gsl_matrix_scale(m_, d);
407    return *this;
408  }
409
410
411  bool isnull(const matrix& m)
412  {
413    return gsl_matrix_isnull(m.gsl_matrix_p());
414  }
415
416
417  double max(const matrix& m)
418  {
419    return gsl_matrix_max(m.gsl_matrix_p());
420  }
421
422
423  double min(const matrix& m)
424  {
425    return gsl_matrix_min(m.gsl_matrix_p());
426  }
427
428
429  void minmax_index(const matrix& m,
430                    std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
431  {
432    gsl_matrix_minmax_index(m.gsl_matrix_p(), &min.first, &min.second,
433                            &max.first, &max.second);
434  }
435
436
437  bool nan(const matrix& templat, matrix& flag)
438  {
439    size_t rows=templat.rows();
440    size_t columns=templat.columns();
441    if (rows!=flag.rows() && columns!=flag.columns())
442      flag=matrix(rows,columns,1.0);
443    else
444      flag.set_all(1.0);
445    bool nan=false;
446    for (size_t i=0; i<rows; i++)
447      for (size_t j=0; j<columns; j++) 
448        if (std::isnan(templat(i,j))) {
449          flag(i,j)=0;
450          nan=true;
451        }
452    return nan;
453  }
454
455
456  void swap(matrix& a, matrix& b)
457  {
458    int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
459    if (status)
460      throw utility::GSL_error(std::string("swap(matrix&,matrix&)",status));
461  }
462
463
464  std::ostream& operator<<(std::ostream& s, const matrix& m)
465  {
466    s.setf(std::ios::dec);
467    s.precision(12);
468    for(size_t i=0, j=0; i<m.rows(); i++)
469      for (j=0; j<m.columns(); j++) {
470        s << m(i,j);
471        if (j<m.columns()-1)
472          s << s.fill();
473        else if (i<m.rows()-1)
474          s << "\n";
475      }
476    return s;
477  }
478
479
480  vector operator*(const matrix& m, const vector& v)
481  {
482    utility::vector res(m.rows());
483    for (size_t i=0; i<res.size(); ++i)
484      res(i) = vector(m,i) * v;
485    return res;
486  }
487
488
489  vector operator*(const vector& v, const matrix& m)
490  {
491    utility::vector res(m.columns());
492    for (size_t i=0; i<res.size(); ++i)
493      res(i) = v * vector(m,i,false);
494    return res;
495  }
496
497}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.