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

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

fixes #42

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.0 KB
Line 
1// $Id: matrix.cc 732 2007-01-06 16:11:22Z 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  const vector matrix::operator*(const vector&) const
362  {
363    std::cerr << "Not implemented:" << std::endl
364              << "   const vector matrix::operator*(const vector&) const"
365              << std::endl;
366    return vector(0);
367  }
368
369
370  bool matrix::operator==(const matrix& other) const
371  {
372    return equal(other);
373  }
374
375
376  bool matrix::operator!=(const matrix& other) const
377  {
378    return !equal(other);
379  }
380
381
382  const matrix& matrix::operator=( const matrix& other )
383  {
384    if ( this != &other ) {
385      if (view_)
386        delete view_;
387      else if (m_)
388        gsl_matrix_free(m_);
389      m_ = other.create_gsl_matrix_copy();
390    }
391    return *this;
392  }
393
394
395  const matrix& matrix::operator+=(const matrix& m)
396  {
397    add(m);
398    return *this;
399  }
400
401
402  const matrix& matrix::operator+=(const double d)
403  {
404    add_constant(d);
405    return *this;
406  }
407
408
409  const matrix& matrix::operator-=(const matrix& m)
410  {
411    sub(m);
412    return *this;
413  }
414
415
416  const matrix& matrix::operator*=(const matrix& other)
417  {
418    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
419    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
420    gsl_matrix_free(m_);
421    m_=result;
422    return *this;
423  }
424
425
426  const matrix& matrix::operator*=(const double d)
427  {
428    scale(d);
429    return *this;
430  }
431
432
433  std::ostream& operator<<(std::ostream& s, const matrix& m)
434  {
435    s.setf(std::ios::dec);
436    s.precision(12);
437    for(size_t i=0, j=0; i<m.rows(); i++)
438      for (j=0; j<m.columns(); j++) {
439        s << m(i,j);
440        if (j<m.columns()-1)
441          s << s.fill();
442        else if (i<m.rows()-1)
443          s << "\n";
444      }
445    return s;
446  }
447
448}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.