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

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

Fixed exception message.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.0 KB
Line 
1// $Id: matrix.cc 780 2007-03-05 19:30:49Z 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 (columns()!=other.columns() || rows()!=other.rows())
192      return false;
193    for (size_t i=0; i<rows(); i++)
194      for (size_t j=0; j<columns(); j++)
195        // The two last condition checks are needed for NaN detection
196        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
197            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
198          return false;
199    return true;
200  }
201
202
203  const gsl_matrix* matrix::gsl_matrix_p(void) const
204  {
205    return m_;
206  }
207
208
209  gsl_matrix* matrix::gsl_matrix_p(void)
210  {
211    return m_;
212  }
213
214
215  bool matrix::isview(void) const
216  {
217    return view_;
218  }
219
220
221  void matrix::mul_elements(const matrix& b)
222  {
223    int status=gsl_matrix_mul_elements(m_, b.m_);
224    if (status)
225      throw utility::GSL_error(std::string("matrix::mul_elements",status));
226   
227  }
228
229
230  size_t matrix::rows(void) const
231  {
232    if (!m_)
233      return 0;
234    return m_->size1;
235  }
236
237
238  void matrix::set(const matrix& mat)
239  {
240    if (gsl_matrix_memcpy(m_, mat.m_))
241      throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mis-match");
242  }
243
244
245  void matrix::set_all(const double value)
246  {
247    gsl_matrix_set_all(m_, value);
248  }
249
250
251  void matrix::set_column(const size_t column, const vector& vec)
252  {
253    int status=gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
254    if (status)
255      throw utility::GSL_error(std::string("matrix::set_column",status));
256  }
257
258
259  void matrix::set_row(const size_t row, const vector& vec)
260  {
261    int status=gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
262    if (status)
263      throw utility::GSL_error(std::string("matrix::set_row",status));
264  }
265
266
267  void matrix::swap_columns(const size_t i, const size_t j)
268  {
269    int status=gsl_matrix_swap_columns(m_, i, j);
270    if (status)
271      throw utility::GSL_error(std::string("matrix::swap_columns",status));
272  }
273
274
275  void matrix::swap_rowcol(const size_t i, const size_t j)
276  {
277    int status=gsl_matrix_swap_rowcol(m_, i, j);
278    if (status)
279      throw utility::GSL_error(std::string("matrix::swap_rowcol",status));
280  }
281
282
283  void matrix::swap_rows(const size_t i, const size_t j)
284  {
285    int status=gsl_matrix_swap_rows(m_, i, j);
286    if (status)
287      throw utility::GSL_error(std::string("matrix::swap_rows",status));
288  }
289
290
291  void matrix::transpose(void)
292  {
293    if (columns()==rows())
294      gsl_matrix_transpose(m_); // this never fails
295    else {
296      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
297      if (!transposed)
298        throw utility::GSL_error("matrix::transpose failed to allocate memory");
299      // next line never fails if allocation above succeeded.
300      gsl_matrix_transpose_memcpy(transposed,m_);
301      gsl_matrix_free(m_);
302      m_=transposed;
303      if (blas_result_) {
304        gsl_matrix_free(blas_result_);
305        blas_result_=NULL;
306      }
307    }
308  }
309
310
311  double& matrix::operator()(size_t row, size_t column)
312  {
313    double* d=gsl_matrix_ptr(m_, row, column);
314    if (!d)
315      throw utility::GSL_error("matrix::operator()",GSL_EINVAL);
316    return *d;
317  }
318
319
320  const double& matrix::operator()(size_t row, size_t column) const
321  {
322    const double* d=gsl_matrix_const_ptr(m_, row, column);
323    if (!d)
324      throw utility::GSL_error("matrix::operator()",GSL_EINVAL);
325    return *d;
326  }
327
328
329  bool matrix::operator==(const matrix& other) const
330  {
331    return equal(other);
332  }
333
334
335  bool matrix::operator!=(const matrix& other) const
336  {
337    return !equal(other);
338  }
339
340
341  const matrix& matrix::operator=( const matrix& other )
342  {
343    if ( this != &other ) {
344      if (view_)
345        delete view_;
346      else if (m_)
347        gsl_matrix_free(m_);
348      m_ = other.create_gsl_matrix_copy();
349    }
350    return *this;
351  }
352
353
354  const matrix& matrix::operator+=(const matrix& m)
355  {
356    int status=gsl_matrix_add(m_, m.m_);
357    if (status)
358      throw utility::GSL_error(std::string("matrix::operator+=", status));
359    return *this;
360  }
361
362
363  const matrix& matrix::operator+=(const double d)
364  {
365    gsl_matrix_add_constant(m_, d);
366    return *this;
367  }
368
369
370  const matrix& matrix::operator-=(const matrix& m)
371  {
372    int status=gsl_matrix_sub(m_, m.m_);
373    if (status)
374      throw utility::GSL_error(std::string("matrix::operator-=", status));
375    return *this;
376  }
377
378
379  const matrix& matrix::operator-=(const double d)
380  {
381    gsl_matrix_add_constant(m_, -d);
382    return *this;
383  }
384
385
386  const matrix& matrix::operator*=(const matrix& other)
387  {
388    if (!blas_result_) {
389      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
390      if (!blas_result_)
391        throw utility::GSL_error("matrix::operator*= failed to allocate memory");
392    }
393    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0,
394                   blas_result_);
395    gsl_matrix* tmp=m_;
396    m_=blas_result_;
397    blas_result_=tmp;
398    return *this;
399  }
400
401
402  const matrix& matrix::operator*=(const double d)
403  {
404    gsl_matrix_scale(m_, d);
405    return *this;
406  }
407
408
409  bool isnull(const matrix& m)
410  {
411    return gsl_matrix_isnull(m.gsl_matrix_p());
412  }
413
414
415  double max(const matrix& m)
416  {
417    return gsl_matrix_max(m.gsl_matrix_p());
418  }
419
420
421  double min(const matrix& m)
422  {
423    return gsl_matrix_min(m.gsl_matrix_p());
424  }
425
426
427  void minmax_index(const matrix& m,
428                    std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
429  {
430    gsl_matrix_minmax_index(m.gsl_matrix_p(), &min.first, &min.second,
431                            &max.first, &max.second);
432  }
433
434
435  bool nan(const matrix& templat, matrix& flag)
436  {
437    size_t rows=templat.rows();
438    size_t columns=templat.columns();
439    if (rows!=flag.rows() && columns!=flag.columns())
440      flag=matrix(rows,columns,1.0);
441    else
442      flag.set_all(1.0);
443    bool nan=false;
444    for (size_t i=0; i<rows; i++)
445      for (size_t j=0; j<columns; j++) 
446        if (std::isnan(templat(i,j))) {
447          flag(i,j)=0;
448          nan=true;
449        }
450    return nan;
451  }
452
453
454  void swap(matrix& a, matrix& b)
455  {
456    int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
457    if (status)
458      throw utility::GSL_error(std::string("swap(matrix&,matrix&)",status));
459  }
460
461
462  std::ostream& operator<<(std::ostream& s, const matrix& m)
463  {
464    s.setf(std::ios::dec);
465    s.precision(12);
466    for(size_t i=0, j=0; i<m.rows(); i++)
467      for (j=0; j<m.columns(); j++) {
468        s << m(i,j);
469        if (j<m.columns()-1)
470          s << s.fill();
471        else if (i<m.rows()-1)
472          s << "\n";
473      }
474    return s;
475  }
476
477
478  vector operator*(const matrix& m, const vector& v)
479  {
480    utility::vector res(m.rows());
481    for (size_t i=0; i<res.size(); ++i)
482      res(i) = vector(m,i) * v;
483    return res;
484  }
485
486
487  vector operator*(const vector& v, const matrix& m)
488  {
489    utility::vector res(m.columns());
490    for (size_t i=0; i<res.size(); ++i)
491      res(i) = v * vector(m,i,false);
492    return res;
493  }
494
495}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.