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

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

fixes #306 matrix iterators

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