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

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

Addresses #170.

  • 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 716 2006-12-25 00:53:13Z 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
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 "stl_utility.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    : m_(NULL), view_(NULL)
44  {
45  }
46
47
48  matrix::matrix(const size_t& r, const size_t& c, double init_value)
49    : view_(NULL)
50  {
51    m_ = gsl_matrix_alloc(r,c);
52    set_all(init_value);
53  }
54
55
56  matrix::matrix(const matrix& o)
57    : view_(NULL)
58  {
59    m_ = o.create_gsl_matrix_copy();
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  {
66    // Jari, exception handling needed here. Failure in setting up a
67    // proper gsl_matrix_view is communicated by NULL pointer in the
68    // view structure (cf. GSL manual). How about GSL error state?
69    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
70                                                     offset_row,offset_column,
71                                                     n_row,n_column));
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    : view_(NULL)
80  {
81    // Markus to Jari, somewhere we should check that quiet_NaNs are supported
82    // std::numeric_limits<double>::has_quiet_NaN has to be true.
83    // Also in vector
84
85    // read the data file and store in stl vectors (dynamically
86    // expandable)
87    std::vector<std::vector<double> > data_matrix;
88    u_int nof_columns=0;
89    u_int nof_rows = 0;
90    std::string line;
91    while(getline(is, line, '\n')){
92      // Ignoring empty lines
93      if (!line.size()) {
94        continue;
95      }
96      nof_rows++;
97      std::vector<double> v;
98      std::string element;
99      std::stringstream ss(line);
100     
101      bool ok=true;
102      while(ok) {
103        if(sep=='\0')
104          ok=(ss>>element);
105        else
106          ok=getline(ss, element, sep);
107        if(!ok)
108          break;
109       
110        if(utility::is_double(element)) {
111          v.push_back(atof(element.c_str()));
112        }
113        else if (!element.size() || utility::is_nan(element)) {
114          v.push_back(std::numeric_limits<double>::quiet_NaN());
115        }
116        else {
117          // Jari, this should be communicated with as an exception.
118          std::cerr << "Warning: '" << element
119                    << "' is not accepted as a matrix element." << std::endl;
120        }
121      }           
122      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
123          v.push_back(std::numeric_limits<double>::quiet_NaN());
124      if (!nof_columns)
125        nof_columns=v.size();
126      else if (v.size()!=nof_columns) {
127        std::ostringstream s;
128        s << "matrix::matrix(std::istream&, char) data file error: "
129          << "line " << nof_rows << " has " << v.size()
130          << " columns; expected " << nof_columns << " columns.";
131        throw utility::IO_error(s.str());
132      }
133      data_matrix.push_back(v);
134    }
135
136    // manipulate the state of the stream to be good
137    is.clear(std::ios::goodbit);
138    // convert the data to a gsl matrix
139    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
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    if (view_)
149      delete view_;
150    else if (m_)
151      gsl_matrix_free(m_);
152    m_=NULL;
153  }
154
155
156  int matrix::add(const matrix& b)
157  {
158    return gsl_matrix_add(m_, b.m_);
159  }
160
161
162  int matrix::add_constant(const double d)
163  {
164    return gsl_matrix_add_constant(m_, d);
165  }
166
167
168  size_t matrix::columns(void) const
169  {
170    return m_->size2;
171  }
172
173
174  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
175  {
176    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
177    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
178    return m;
179  }
180
181
182  int matrix::div_elements(const matrix& b)
183  {
184    return gsl_matrix_div_elements(m_, b.m_);
185  }
186
187
188  bool matrix::equal(const matrix& other, const double d) const
189  {
190    if (columns()!=other.columns() || rows()!=other.rows())
191      return false;
192    for (size_t i=0; i<rows(); i++)
193      for (size_t j=0; j<columns(); j++)
194        // The two last condition checks are needed for NaN detection
195        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
196            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
197          return false;
198    return true;
199  }
200
201
202  const gsl_matrix* matrix::gsl_matrix_p(void) const
203  {
204    return m_;
205  }
206
207
208  gsl_matrix* matrix::gsl_matrix_p(void)
209  {
210    return m_;
211  }
212
213
214  bool matrix::isnull(void) const
215  {
216    return gsl_matrix_isnull(m_);
217  }
218
219
220  bool matrix::isview(void) const
221  {
222    return view_;
223  }
224
225
226  double matrix::max(void) const
227  {
228    return gsl_matrix_max(m_);
229  }
230
231
232  double matrix::min(void) const
233  {
234    return gsl_matrix_min(m_);
235  }
236
237
238  void matrix::minmax_index(std::pair<size_t,size_t>& min,
239                            std::pair<size_t,size_t>& max) const
240  {
241    gsl_matrix_minmax_index(m_, &min.first, &min.second, &max.first,
242                            &max.second);
243  }
244
245
246  bool matrix::nan(matrix &m) const
247  {
248    m=matrix(rows(),columns(),1.0);
249    bool nan=false;
250    for (size_t i=0; i<rows(); i++)
251      for (size_t j=0; j<columns(); j++) 
252        if (std::isnan(operator()(i,j))) {
253          m(i,j)=0;
254          nan=true;
255        }
256    return nan;
257  }
258
259
260  int matrix::mul_elements(const matrix& b)
261  {
262    return gsl_matrix_mul_elements(m_, b.m_);
263  }
264
265
266  size_t matrix::rows(void) const
267  {
268    return m_->size1;
269  }
270
271
272  int matrix::scale(const double d)
273  {
274    return gsl_matrix_scale(m_, d);
275  }
276
277
278  int matrix::set(const matrix& mat)
279  {
280    return gsl_matrix_memcpy(m_, mat.m_);
281  }
282
283
284  void matrix::set_all(const double value)
285  {
286    gsl_matrix_set_all(m_, value);
287  }
288
289
290  int matrix::set_column(const size_t column, const vector& vec)
291  {
292    return gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
293  }
294
295
296  int matrix::set_row(const size_t row, const vector& vec)
297  {
298    return gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
299  }
300
301
302  int matrix::sub(const matrix& b)
303  {
304    return gsl_matrix_sub(m_, b.m_);
305  }
306
307
308  int matrix::swap(matrix& other)
309  {
310    return gsl_matrix_swap(m_, other.m_);
311  }
312
313
314  int matrix::swap_columns(const size_t i, const size_t j)
315  {
316    return gsl_matrix_swap_columns(m_, i, j);
317  }
318
319
320  int matrix::swap_rowcol(const size_t i, const size_t j)
321  {
322    return gsl_matrix_swap_rowcol(m_, i, j);
323  }
324
325
326  int matrix::swap_rows(const size_t i, const size_t j)
327  {
328    return gsl_matrix_swap_rows(m_, i, j);
329  }
330
331
332  // Jari, checkout GSL transpose support in GSL manual 8.4.9
333  void matrix::transpose(void)
334  {
335    if (columns()==rows())
336      gsl_matrix_transpose(m_);
337    else {
338      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
339      gsl_matrix_transpose_memcpy(transposed,m_);
340      gsl_matrix_free(m_);
341      m_=transposed;
342    }
343  }
344
345
346  double& matrix::operator()(size_t row, size_t column)
347  {
348    return (*gsl_matrix_ptr(m_, row, column));
349  }
350
351
352  const double& matrix::operator()(size_t row, size_t column) const
353  {
354    return (*gsl_matrix_const_ptr(m_, row, column));
355  }
356
357
358  const vector matrix::operator*(const vector&) const
359  {
360    std::cerr << "Not implemented:" << std::endl
361              << "   const vector matrix::operator*(const vector&) const"
362              << std::endl;
363    return vector(0);
364  }
365
366
367  bool matrix::operator==(const matrix& other) const
368  {
369    return equal(other);
370  }
371
372
373  bool matrix::operator!=(const matrix& other) const
374  {
375    return !equal(other);
376  }
377
378
379  const matrix& matrix::operator=( const matrix& other )
380  {
381    if ( this != &other ) {
382      if (view_)
383        delete view_;
384      else if (m_)
385        gsl_matrix_free(m_);
386      m_ = other.create_gsl_matrix_copy();
387    }
388    return *this;
389  }
390
391
392  const matrix& matrix::operator+=(const matrix& m)
393  {
394    add(m);
395    return *this;
396  }
397
398
399  const matrix& matrix::operator+=(const double d)
400  {
401    add_constant(d);
402    return *this;
403  }
404
405
406  const matrix& matrix::operator-=(const matrix& m)
407  {
408    sub(m);
409    return *this;
410  }
411
412
413  const matrix& matrix::operator*=(const matrix& other)
414  {
415    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
416    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
417    gsl_matrix_free(m_);
418    m_=result;
419    return *this;
420  }
421
422
423  const matrix& matrix::operator*=(const double d)
424  {
425    scale(d);
426    return *this;
427  }
428
429
430  std::ostream& operator<<(std::ostream& s, const matrix& m)
431  {
432    s.setf(std::ios::dec);
433    s.precision(12);
434    for(size_t i=0, j=0; i<m.rows(); i++)
435      for (j=0; j<m.columns(); j++) {
436        s << m(i,j);
437        if (j<m.columns()-1)
438          s << s.fill();
439        else if (i<m.rows()-1)
440          s << "\n";
441      }
442    return s;
443  }
444
445}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.