source: trunk/yat/utility/Matrix.cc @ 1147

Last change on this file since 1147 was 1147, checked in by Peter, 14 years ago

fixes #294

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.2 KB
Line 
1// $Id: Matrix.cc 1147 2008-02-25 21:24:40Z 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, Markus Ringnér, Peter Johansson
7  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://trac.thep.lu.se/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 "VectorBase.h"
30#include "VectorConstView.h"
31#include "VectorView.h"
32#include "utility.h"
33
34#include <cassert>
35#include <cmath>
36#include <sstream>
37#include <vector>
38
39#include <gsl/gsl_blas.h>
40
41namespace theplu {
42namespace yat {
43namespace utility {
44
45
46  Matrix::Matrix(void)
47    : blas_result_(NULL), m_(NULL)
48  {
49  }
50
51
52  Matrix::Matrix(const size_t& r, const size_t& c, double init_value)
53    : blas_result_(NULL), m_(gsl_matrix_alloc(r,c))
54  {
55    if (!m_)
56      throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
57    all(init_value);
58  }
59
60
61  Matrix::Matrix(const Matrix& o)
62    : blas_result_(NULL), m_(o.create_gsl_matrix_copy())
63  {
64  }
65
66
67  // Constructor that gets data from istream
68  Matrix::Matrix(std::istream& is, char sep) 
69    throw (utility::IO_error,std::exception)
70    : blas_result_(NULL)
71  {
72    if (!is.good())
73      throw utility::IO_error("Matrix: istream is not good");
74
75    // read the data file and store in stl vectors (dynamically
76    // expandable)
77    std::vector<std::vector<double> > data_matrix;
78    u_int nof_columns=0;
79    u_int nof_rows = 0;
80    std::string line;
81    while(getline(is, line, '\n')){
82      // Ignoring empty lines
83      if (!line.size()) {
84        continue;
85      }
86      nof_rows++;
87      std::vector<double> v;
88      std::string element;
89      std::stringstream ss(line);
90     
91      bool ok=true;
92      while(ok) {
93        if(sep=='\0')
94          ok=(ss>>element);
95        else
96          ok=getline(ss, element, sep);
97        if(!ok)
98          break;
99       
100        if(utility::is_double(element)) {
101          v.push_back(atof(element.c_str()));
102        }
103        else if (!element.size() || utility::is_nan(element)) {
104          v.push_back(std::numeric_limits<double>::quiet_NaN());
105        }
106        else {
107          std::stringstream ss("Warning: '");
108          ss << element << "' is not accepted as a matrix element.";
109          throw IO_error(ss.str());
110        }
111      }           
112      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
113          v.push_back(std::numeric_limits<double>::quiet_NaN());
114      if (!nof_columns)
115        nof_columns=v.size();
116      else if (v.size()!=nof_columns) {
117        std::ostringstream s;
118        s << "Matrix::Matrix(std::istream&, char) data file error: "
119          << "line " << nof_rows << " has " << v.size()
120          << " columns; expected " << nof_columns << " columns.";
121        throw utility::IO_error(s.str());
122      }
123      data_matrix.push_back(v);
124    }
125
126    // manipulate the state of the stream to be good
127    is.clear(std::ios::goodbit);
128   
129    // if stream was empty, create nothing
130    if (!nof_columns || !nof_rows)
131      return;
132
133    // convert the data to a gsl matrix
134    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
135    if (!m_)
136      throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
137
138    // if gsl error handler disabled, out of bounds index will not
139    // abort the program.
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    delete_allocated_memory();
149    if (blas_result_)
150      gsl_matrix_free(blas_result_);
151    blas_result_=NULL;
152  }
153
154
155  void Matrix::all(const double value)
156  {
157    assert(m_);
158    gsl_matrix_set_all(m_, value);
159  }
160
161
162  Matrix::iterator Matrix::begin(void)
163  {
164    return iterator(&(*this)(0,0), 1);
165  }
166
167
168  Matrix::const_iterator Matrix::begin(void) const
169  {
170    return const_iterator(&(*this)(0,0), 1);
171  }
172
173
174  Matrix::column_iterator Matrix::begin_column(size_t i)
175  {
176    return iterator(&(*this)(0,i), this->columns());
177  }
178
179
180  Matrix::const_column_iterator Matrix::begin_column(size_t i) const
181  {
182    return const_iterator(&(*this)(0,i), this->columns());
183  }
184
185
186  Matrix::row_iterator Matrix::begin_row(size_t i)
187  {
188    return iterator(&(*this)(i,0), 1);
189  }
190
191
192  Matrix::const_row_iterator Matrix::begin_row(size_t i) const
193  {
194    return const_iterator(&(*this)(i,0), 1);
195  }
196
197
198  VectorView Matrix::column_view(size_t col)
199  {
200    VectorView res(*this, col, false);
201    return res;
202  }
203
204
205  const VectorConstView Matrix::column_const_view(size_t col) const
206  {
207    return VectorConstView(*this, col, false);
208  }
209
210
211  size_t Matrix::columns(void) const
212  {
213    return (m_ ? m_->size2 : 0);
214  }
215
216
217  gsl_matrix* Matrix::create_gsl_matrix_copy(void) const
218  {
219    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
220    if (!m)
221      throw utility::GSL_error("Matrix::create_gsl_matrix_copy failed to allocate memory");
222    if (gsl_matrix_memcpy(m,m_))
223      throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
224    return m;
225  }
226
227
228  void Matrix::delete_allocated_memory(void)
229  {
230    if (m_)
231      gsl_matrix_free(m_);
232    m_=NULL;
233  }
234
235
236  void Matrix::div(const Matrix& other)
237  {
238    assert(m_);
239    int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
240    if (status)
241      throw utility::GSL_error(std::string("matrix::div_elements",status));
242  }
243
244
245  Matrix::iterator Matrix::end(void)
246  {
247    return iterator(&(*this)(0,0)+rows()*columns(), 1);
248  }
249
250
251  Matrix::const_iterator Matrix::end(void) const
252  {
253    return const_iterator(&(*this)(0,0)+rows()*columns(), 1);
254  }
255
256
257  Matrix::column_iterator Matrix::end_column(size_t i)
258  {
259    return column_iterator(&(*this)(0,i)+rows()*columns(), this->columns());
260  }
261
262
263  Matrix::const_column_iterator Matrix::end_column(size_t i) const
264  {
265    return const_column_iterator(&(*this)(0,i)+rows()*columns(),this->columns());
266  }
267
268
269  Matrix::row_iterator Matrix::end_row(size_t i)
270  {
271    return row_iterator(&(*this)(i,0)+columns(), 1);
272  }
273
274
275  Matrix::const_row_iterator Matrix::end_row(size_t i) const
276  {
277    return const_row_iterator(&(*this)(i,0)+columns(), 1);
278  }
279
280
281  bool Matrix::equal(const Matrix& other, const double d) const
282  {
283    if (this==&other)
284      return true;
285    if (columns()!=other.columns() || rows()!=other.rows())
286      return false;
287    for (size_t i=0; i<rows(); i++)
288      for (size_t j=0; j<columns(); j++)
289        // The two last condition checks are needed for NaN detection
290        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
291            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
292          return false;
293    return true;
294  }
295
296
297  const gsl_matrix* Matrix::gsl_matrix_p(void) const
298  {
299    return m_;
300  }
301
302
303  gsl_matrix* Matrix::gsl_matrix_p(void)
304  {
305    return m_;
306  }
307
308
309  void Matrix::mul(const Matrix& other)
310  {
311    assert(m_);
312    int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());
313    if (status)
314      throw utility::GSL_error(std::string("Matrix::mul_elements",status));
315  }
316
317
318  void Matrix::resize(size_t r, size_t c, double init_value)
319  {
320    delete_allocated_memory();
321
322    m_ = gsl_matrix_alloc(r,c);
323    if (!m_)
324      throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
325    all(init_value);
326
327    // no need to delete blas_result_ if the number of rows fit, it
328    // may be useful later.
329    if (blas_result_ && (blas_result_->size1!=rows())) {
330      gsl_matrix_free(blas_result_);
331      blas_result_=NULL;
332    }
333  }
334
335
336  size_t Matrix::rows(void) const
337  {
338    return (m_ ? m_->size1 : 0);
339  }
340
341
342  const VectorConstView Matrix::row_const_view(size_t col) const
343  {
344    return VectorConstView(*this, col, true);
345  }
346
347
348  VectorView Matrix::row_view(size_t row)
349  {
350    VectorView res(*this, row, true);
351    return res;
352  }
353
354
355  void Matrix::swap_columns(const size_t i, const size_t j)
356  {
357    assert(m_);
358    int status=gsl_matrix_swap_columns(m_, i, j);
359    if (status)
360      throw utility::GSL_error(std::string("Matrix::swap_columns",status));
361  }
362
363
364  void Matrix::swap_rowcol(const size_t i, const size_t j)
365  {
366    assert(m_);
367    int status=gsl_matrix_swap_rowcol(m_, i, j);
368    if (status)
369      throw utility::GSL_error(std::string("Matrix::swap_rowcol",status));
370  }
371
372
373  void Matrix::swap_rows(const size_t i, const size_t j)
374  {
375    assert(m_);
376    int status=gsl_matrix_swap_rows(m_, i, j);
377    if (status)
378      throw utility::GSL_error(std::string("Matrix::swap_rows",status));
379  }
380
381
382  void Matrix::transpose(void)
383  {
384    assert(m_);
385    if (columns()==rows())
386      gsl_matrix_transpose(m_); // this never fails
387    else {
388      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
389      if (!transposed)
390        throw utility::GSL_error("Matrix::transpose failed to allocate memory");
391      // next line never fails if allocation above succeeded.
392      gsl_matrix_transpose_memcpy(transposed,m_);
393      gsl_matrix_free(m_);
394      m_ = transposed;
395      if (blas_result_) {
396        gsl_matrix_free(blas_result_);
397        blas_result_=NULL;
398      }
399    }
400  }
401
402
403  double& Matrix::operator()(size_t row, size_t column)
404  {
405    assert(m_);
406    assert(row<rows());
407    assert(column<columns());
408    double* d=gsl_matrix_ptr(m_, row, column);
409    if (!d)
410      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
411    return *d;
412  }
413
414
415  const double& Matrix::operator()(size_t row, size_t column) const
416  {
417    assert(row<rows());
418    assert(column<columns());
419    const double* d=gsl_matrix_const_ptr(m_, row, column);
420    if (!d)
421      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
422    return *d;
423  }
424
425
426  bool Matrix::operator==(const Matrix& other) const
427  {
428    return equal(other);
429  }
430
431
432  bool Matrix::operator!=(const Matrix& other) const
433  {
434    return !equal(other);
435  }
436
437
438  const Matrix& Matrix::operator=( const Matrix& other )
439  {
440    assert(other.m_);
441    if (this!=&other) {
442      if ( !m_ || (other.m_->size1!=m_->size1) || (other.m_->size2!=m_->size2) )
443        resize(other.m_->size1,other.m_->size2);
444      if (gsl_matrix_memcpy(m_, other.gsl_matrix_p()))
445        throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
446    }
447    return *this;
448  }
449
450
451  const Matrix& Matrix::operator+=(const Matrix& other)
452  {
453    assert(m_);
454    int status=gsl_matrix_add(m_, other.m_);
455    if (status)
456      throw utility::GSL_error(std::string("Matrix::operator+=", status));
457    return *this;
458  }
459
460
461  const Matrix& Matrix::operator+=(const double d)
462  {
463    assert(m_);
464    gsl_matrix_add_constant(m_, d);
465    return *this;
466  }
467
468
469  const Matrix& Matrix::operator-=(const Matrix& other)
470  {
471    assert(m_);
472    int status=gsl_matrix_sub(m_, other.m_);
473    if (status)
474      throw utility::GSL_error(std::string("Matrix::operator-=", status));
475    return *this;
476  }
477
478
479  const Matrix& Matrix::operator-=(const double d)
480  {
481    assert(m_);
482    gsl_matrix_add_constant(m_, -d);
483    return *this;
484  }
485
486
487  const Matrix& Matrix::operator*=(const Matrix& other)
488  {
489    assert(m_);
490    if ( blas_result_ && ((blas_result_->size1!=rows()) ||
491                          (blas_result_->size2!=other.columns())) ) {
492      gsl_matrix_free(blas_result_);
493      blas_result_=NULL;
494    }
495    if (!blas_result_) {
496      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
497      if (!blas_result_)
498        throw utility::GSL_error("Matrix::operator*= failed to allocate memory");
499    }
500    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, blas_result_);
501    gsl_matrix* tmp=m_;
502    m_ = blas_result_;
503    blas_result_=tmp;
504    return *this;
505  }
506
507
508  const Matrix& Matrix::operator*=(const double d)
509  {
510    assert(m_);
511    gsl_matrix_scale(m_, d);
512    return *this;
513  }
514
515
516  bool isnull(const Matrix& other)
517  {
518    return gsl_matrix_isnull(other.gsl_matrix_p());
519  }
520
521
522  double max(const Matrix& other)
523  {
524    return gsl_matrix_max(other.gsl_matrix_p());
525  }
526
527
528  double min(const Matrix& other)
529  {
530    return gsl_matrix_min(other.gsl_matrix_p());
531  }
532
533
534  void minmax_index(const Matrix& other,
535                    std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
536  {
537    gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,
538                            &max.first, &max.second);
539  }
540
541
542  bool nan(const Matrix& templat, Matrix& flag)
543  {
544    size_t rows=templat.rows();
545    size_t columns=templat.columns();
546    if (rows!=flag.rows() && columns!=flag.columns())
547      flag.resize(rows,columns,1.0);
548    else
549      flag.all(1.0);
550    bool nan=false;
551    for (size_t i=0; i<rows; i++)
552      for (size_t j=0; j<columns; j++) 
553        if (std::isnan(templat(i,j))) {
554          flag(i,j)=0;
555          nan=true;
556        }
557    return nan;
558  }
559
560
561  void swap(Matrix& a, Matrix& b)
562  {
563    assert(a.gsl_matrix_p()); assert(b.gsl_matrix_p());
564    int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
565    if (status)
566      throw utility::GSL_error(std::string("swap(Matrix&,Matrix&)",status));
567  }
568
569
570  std::ostream& operator<<(std::ostream& s, const Matrix& m)
571  {
572    s.setf(std::ios::dec);
573    s.precision(12);
574    for(size_t i=0, j=0; i<m.rows(); i++)
575      for (j=0; j<m.columns(); j++) {
576        s << m(i,j);
577        if (j<m.columns()-1)
578          s << s.fill();
579        else if (i<m.rows()-1)
580          s << "\n";
581      }
582    return s;
583  }
584
585
586  Vector operator*(const Matrix& m, const VectorBase& v)
587  {
588    utility::Vector res(m.rows());
589    for (size_t i=0; i<res.size(); ++i)
590      res(i) = VectorConstView(m,i) * v;
591    return res;
592  }
593
594
595  Vector operator*(const VectorBase& v, const Matrix& m)
596  {
597    utility::Vector res(m.columns());
598    for (size_t i=0; i<res.size(); ++i)
599      res(i) = v * VectorConstView(m,i,false);
600    return res;
601  }
602
603}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.