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

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

removing includes

  • 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 725 2007-01-01 21:17:41Z 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    return m_->size2;
170  }
171
172
173  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
174  {
175    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
176    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
177    return m;
178  }
179
180
181  int matrix::div_elements(const matrix& b)
182  {
183    return gsl_matrix_div_elements(m_, b.m_);
184  }
185
186
187  bool matrix::equal(const matrix& other, const double d) const
188  {
189    if (columns()!=other.columns() || rows()!=other.rows())
190      return false;
191    for (size_t i=0; i<rows(); i++)
192      for (size_t j=0; j<columns(); j++)
193        // The two last condition checks are needed for NaN detection
194        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
195            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
196          return false;
197    return true;
198  }
199
200
201  const gsl_matrix* matrix::gsl_matrix_p(void) const
202  {
203    return m_;
204  }
205
206
207  gsl_matrix* matrix::gsl_matrix_p(void)
208  {
209    return m_;
210  }
211
212
213  bool matrix::isnull(void) const
214  {
215    return gsl_matrix_isnull(m_);
216  }
217
218
219  bool matrix::isview(void) const
220  {
221    return view_;
222  }
223
224
225  double matrix::max(void) const
226  {
227    return gsl_matrix_max(m_);
228  }
229
230
231  double matrix::min(void) const
232  {
233    return gsl_matrix_min(m_);
234  }
235
236
237  void matrix::minmax_index(std::pair<size_t,size_t>& min,
238                            std::pair<size_t,size_t>& max) const
239  {
240    gsl_matrix_minmax_index(m_, &min.first, &min.second, &max.first,
241                            &max.second);
242  }
243
244
245  bool matrix::nan(matrix &m) const
246  {
247    m=matrix(rows(),columns(),1.0);
248    bool nan=false;
249    for (size_t i=0; i<rows(); i++)
250      for (size_t j=0; j<columns(); j++) 
251        if (std::isnan(operator()(i,j))) {
252          m(i,j)=0;
253          nan=true;
254        }
255    return nan;
256  }
257
258
259  int matrix::mul_elements(const matrix& b)
260  {
261    return gsl_matrix_mul_elements(m_, b.m_);
262  }
263
264
265  size_t matrix::rows(void) const
266  {
267    return m_->size1;
268  }
269
270
271  int matrix::scale(const double d)
272  {
273    return gsl_matrix_scale(m_, d);
274  }
275
276
277  int matrix::set(const matrix& mat)
278  {
279    return gsl_matrix_memcpy(m_, mat.m_);
280  }
281
282
283  void matrix::set_all(const double value)
284  {
285    gsl_matrix_set_all(m_, value);
286  }
287
288
289  int matrix::set_column(const size_t column, const vector& vec)
290  {
291    return gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
292  }
293
294
295  int matrix::set_row(const size_t row, const vector& vec)
296  {
297    return gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
298  }
299
300
301  int matrix::sub(const matrix& b)
302  {
303    return gsl_matrix_sub(m_, b.m_);
304  }
305
306
307  int matrix::swap(matrix& other)
308  {
309    return gsl_matrix_swap(m_, other.m_);
310  }
311
312
313  int matrix::swap_columns(const size_t i, const size_t j)
314  {
315    return gsl_matrix_swap_columns(m_, i, j);
316  }
317
318
319  int matrix::swap_rowcol(const size_t i, const size_t j)
320  {
321    return gsl_matrix_swap_rowcol(m_, i, j);
322  }
323
324
325  int matrix::swap_rows(const size_t i, const size_t j)
326  {
327    return gsl_matrix_swap_rows(m_, i, j);
328  }
329
330
331  // Jari, checkout GSL transpose support in GSL manual 8.4.9
332  void matrix::transpose(void)
333  {
334    if (columns()==rows())
335      gsl_matrix_transpose(m_);
336    else {
337      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
338      gsl_matrix_transpose_memcpy(transposed,m_);
339      gsl_matrix_free(m_);
340      m_=transposed;
341    }
342  }
343
344
345  double& matrix::operator()(size_t row, size_t column)
346  {
347    return (*gsl_matrix_ptr(m_, row, column));
348  }
349
350
351  const double& matrix::operator()(size_t row, size_t column) const
352  {
353    return (*gsl_matrix_const_ptr(m_, row, column));
354  }
355
356
357  const vector matrix::operator*(const vector&) const
358  {
359    std::cerr << "Not implemented:" << std::endl
360              << "   const vector matrix::operator*(const vector&) const"
361              << std::endl;
362    return vector(0);
363  }
364
365
366  bool matrix::operator==(const matrix& other) const
367  {
368    return equal(other);
369  }
370
371
372  bool matrix::operator!=(const matrix& other) const
373  {
374    return !equal(other);
375  }
376
377
378  const matrix& matrix::operator=( const matrix& other )
379  {
380    if ( this != &other ) {
381      if (view_)
382        delete view_;
383      else if (m_)
384        gsl_matrix_free(m_);
385      m_ = other.create_gsl_matrix_copy();
386    }
387    return *this;
388  }
389
390
391  const matrix& matrix::operator+=(const matrix& m)
392  {
393    add(m);
394    return *this;
395  }
396
397
398  const matrix& matrix::operator+=(const double d)
399  {
400    add_constant(d);
401    return *this;
402  }
403
404
405  const matrix& matrix::operator-=(const matrix& m)
406  {
407    sub(m);
408    return *this;
409  }
410
411
412  const matrix& matrix::operator*=(const matrix& other)
413  {
414    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
415    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
416    gsl_matrix_free(m_);
417    m_=result;
418    return *this;
419  }
420
421
422  const matrix& matrix::operator*=(const double d)
423  {
424    scale(d);
425    return *this;
426  }
427
428
429  std::ostream& operator<<(std::ostream& s, const matrix& m)
430  {
431    s.setf(std::ios::dec);
432    s.precision(12);
433    for(size_t i=0, j=0; i<m.rows(); i++)
434      for (j=0; j<m.columns(); j++) {
435        s << m(i,j);
436        if (j<m.columns()-1)
437          s << s.fill();
438        else if (i<m.rows()-1)
439          s << "\n";
440      }
441    return s;
442  }
443
444}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.