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

Last change on this file since 734 was 734, checked in by Peter, 15 years ago

added vector*matrix and matrix*vector operators. fixes #18

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.8 KB
Line 
1// $Id: matrix.cc 734 2007-01-06 17:04:54Z 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
8  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
9
10  The yat library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU General Public License as
12  published by the Free Software Foundation; either version 2 of the
13  License, or (at your option) any later version.
14
15  The yat library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  General Public License for more details.
19
20  You should have received a copy of the GNU General Public License
21  along with this program; if not, write to the Free Software
22  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
23  02111-1307, USA.
24*/
25
26#include "matrix.h"
27#include "vector.h"
28#include "utility.h"
29
30#include <cmath>
31#include <sstream>
32#include <vector>
33
34#include <gsl/gsl_blas.h>
35
36namespace theplu {
37namespace yat {
38namespace utility {
39
40
41  matrix::matrix(void)
42    : m_(NULL), view_(NULL)
43  {
44  }
45
46
47  matrix::matrix(const size_t& r, const size_t& c, double init_value)
48    : view_(NULL)
49  {
50    m_ = gsl_matrix_alloc(r,c);
51    set_all(init_value);
52  }
53
54
55  matrix::matrix(const matrix& o)
56    : view_(NULL)
57  {
58    m_ = o.create_gsl_matrix_copy();
59  }
60
61
62  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
63                 size_t n_row, size_t n_column)
64  {
65    // Jari, exception handling needed here. Failure in setting up a
66    // proper gsl_matrix_view is communicated by NULL pointer in the
67    // view structure (cf. GSL manual). How about GSL error state?
68    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
69                                                     offset_row,offset_column,
70                                                     n_row,n_column));
71    m_ = &(view_->matrix);
72  }
73
74
75  // Constructor that gets data from istream
76  matrix::matrix(std::istream& is, char sep) 
77    throw (utility::IO_error,std::exception)
78    : view_(NULL)
79  {
80    // Markus to Jari, somewhere we should check that quiet_NaNs are supported
81    // std::numeric_limits<double>::has_quiet_NaN has to be true.
82    // Also in vector
83
84    // read the data file and store in stl vectors (dynamically
85    // expandable)
86    std::vector<std::vector<double> > data_matrix;
87    u_int nof_columns=0;
88    u_int nof_rows = 0;
89    std::string line;
90    while(getline(is, line, '\n')){
91      // Ignoring empty lines
92      if (!line.size()) {
93        continue;
94      }
95      nof_rows++;
96      std::vector<double> v;
97      std::string element;
98      std::stringstream ss(line);
99     
100      bool ok=true;
101      while(ok) {
102        if(sep=='\0')
103          ok=(ss>>element);
104        else
105          ok=getline(ss, element, sep);
106        if(!ok)
107          break;
108       
109        if(utility::is_double(element)) {
110          v.push_back(atof(element.c_str()));
111        }
112        else if (!element.size() || utility::is_nan(element)) {
113          v.push_back(std::numeric_limits<double>::quiet_NaN());
114        }
115        else {
116          // Jari, this should be communicated with as an exception.
117          std::cerr << "Warning: '" << element
118                    << "' is not accepted as a matrix element." << std::endl;
119        }
120      }           
121      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
122          v.push_back(std::numeric_limits<double>::quiet_NaN());
123      if (!nof_columns)
124        nof_columns=v.size();
125      else if (v.size()!=nof_columns) {
126        std::ostringstream s;
127        s << "matrix::matrix(std::istream&, char) data file error: "
128          << "line " << nof_rows << " has " << v.size()
129          << " columns; expected " << nof_columns << " columns.";
130        throw utility::IO_error(s.str());
131      }
132      data_matrix.push_back(v);
133    }
134
135    // manipulate the state of the stream to be good
136    is.clear(std::ios::goodbit);
137    // convert the data to a gsl matrix
138    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
139    for(u_int i=0;i<nof_rows;i++)
140      for(u_int j=0;j<nof_columns;j++)
141        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
142  }
143
144
145  matrix::~matrix(void)
146  {
147    if (view_)
148      delete view_;
149    else if (m_)
150      gsl_matrix_free(m_);
151    m_=NULL;
152  }
153
154
155  int matrix::add(const matrix& b)
156  {
157    return gsl_matrix_add(m_, b.m_);
158  }
159
160
161  int matrix::add_constant(const double d)
162  {
163    return gsl_matrix_add_constant(m_, d);
164  }
165
166
167  size_t matrix::columns(void) const
168  {
169    if (!m_)
170      return 0;
171    return m_->size2;
172  }
173
174
175  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
176  {
177    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
178    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
179    return m;
180  }
181
182
183  int matrix::div_elements(const matrix& b)
184  {
185    return gsl_matrix_div_elements(m_, b.m_);
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::isnull(void) const
216  {
217    return gsl_matrix_isnull(m_);
218  }
219
220
221  bool matrix::isview(void) const
222  {
223    return view_;
224  }
225
226
227  double matrix::max(void) const
228  {
229    return gsl_matrix_max(m_);
230  }
231
232
233  double matrix::min(void) const
234  {
235    return gsl_matrix_min(m_);
236  }
237
238
239  void matrix::minmax_index(std::pair<size_t,size_t>& min,
240                            std::pair<size_t,size_t>& max) const
241  {
242    gsl_matrix_minmax_index(m_, &min.first, &min.second, &max.first,
243                            &max.second);
244  }
245
246
247  bool matrix::nan(matrix &m) const
248  {
249    m=matrix(rows(),columns(),1.0);
250    bool nan=false;
251    for (size_t i=0; i<rows(); i++)
252      for (size_t j=0; j<columns(); j++) 
253        if (std::isnan(operator()(i,j))) {
254          m(i,j)=0;
255          nan=true;
256        }
257    return nan;
258  }
259
260
261  int matrix::mul_elements(const matrix& b)
262  {
263    return gsl_matrix_mul_elements(m_, b.m_);
264  }
265
266
267  size_t matrix::rows(void) const
268  {
269    if (!m_)
270      return 0;
271    return m_->size1;
272  }
273
274
275  int matrix::scale(const double d)
276  {
277    return gsl_matrix_scale(m_, d);
278  }
279
280
281  int matrix::set(const matrix& mat)
282  {
283    return gsl_matrix_memcpy(m_, mat.m_);
284  }
285
286
287  void matrix::set_all(const double value)
288  {
289    gsl_matrix_set_all(m_, value);
290  }
291
292
293  int matrix::set_column(const size_t column, const vector& vec)
294  {
295    return gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
296  }
297
298
299  int matrix::set_row(const size_t row, const vector& vec)
300  {
301    return gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
302  }
303
304
305  int matrix::sub(const matrix& b)
306  {
307    return gsl_matrix_sub(m_, b.m_);
308  }
309
310
311  int matrix::swap(matrix& other)
312  {
313    return gsl_matrix_swap(m_, other.m_);
314  }
315
316
317  int matrix::swap_columns(const size_t i, const size_t j)
318  {
319    return gsl_matrix_swap_columns(m_, i, j);
320  }
321
322
323  int matrix::swap_rowcol(const size_t i, const size_t j)
324  {
325    return gsl_matrix_swap_rowcol(m_, i, j);
326  }
327
328
329  int matrix::swap_rows(const size_t i, const size_t j)
330  {
331    return gsl_matrix_swap_rows(m_, i, j);
332  }
333
334
335  // Jari, checkout GSL transpose support in GSL manual 8.4.9
336  void matrix::transpose(void)
337  {
338    if (columns()==rows())
339      gsl_matrix_transpose(m_);
340    else {
341      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
342      gsl_matrix_transpose_memcpy(transposed,m_);
343      gsl_matrix_free(m_);
344      m_=transposed;
345    }
346  }
347
348
349  double& matrix::operator()(size_t row, size_t column)
350  {
351    return (*gsl_matrix_ptr(m_, row, column));
352  }
353
354
355  const double& matrix::operator()(size_t row, size_t column) const
356  {
357    return (*gsl_matrix_const_ptr(m_, row, column));
358  }
359
360
361  bool matrix::operator==(const matrix& other) const
362  {
363    return equal(other);
364  }
365
366
367  bool matrix::operator!=(const matrix& other) const
368  {
369    return !equal(other);
370  }
371
372
373  const matrix& matrix::operator=( const matrix& other )
374  {
375    if ( this != &other ) {
376      if (view_)
377        delete view_;
378      else if (m_)
379        gsl_matrix_free(m_);
380      m_ = other.create_gsl_matrix_copy();
381    }
382    return *this;
383  }
384
385
386  const matrix& matrix::operator+=(const matrix& m)
387  {
388    add(m);
389    return *this;
390  }
391
392
393  const matrix& matrix::operator+=(const double d)
394  {
395    add_constant(d);
396    return *this;
397  }
398
399
400  const matrix& matrix::operator-=(const matrix& m)
401  {
402    sub(m);
403    return *this;
404  }
405
406
407  const matrix& matrix::operator*=(const matrix& other)
408  {
409    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
410    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
411    gsl_matrix_free(m_);
412    m_=result;
413    return *this;
414  }
415
416
417  const matrix& matrix::operator*=(const double d)
418  {
419    scale(d);
420    return *this;
421  }
422
423
424  std::ostream& operator<<(std::ostream& s, const matrix& m)
425  {
426    s.setf(std::ios::dec);
427    s.precision(12);
428    for(size_t i=0, j=0; i<m.rows(); i++)
429      for (j=0; j<m.columns(); j++) {
430        s << m(i,j);
431        if (j<m.columns()-1)
432          s << s.fill();
433        else if (i<m.rows()-1)
434          s << "\n";
435      }
436    return s;
437  }
438
439}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.