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

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

Fixes #299. Memory leak in matrix was found and removed.

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