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

Last change on this file since 1789 was 1789, checked in by Peter, 12 years ago

Speeding up copying from Matrix to MatrixWeighted? using
BinaryWeight?. Moving BinaryWeight? from stl_utility.h to
utility.h. Improved docs in MatrixWeighted?.

  • 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 1789 2009-02-10 16:11:04Z 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  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://dev.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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include "Matrix.h"
26#include "Vector.h"
27#include "VectorBase.h"
28#include "VectorConstView.h"
29#include "VectorView.h"
30#include "utility.h"
31
32#include <cassert>
33#include <climits>
34#include <cmath>
35#include <sstream>
36#include <vector>
37
38#include <gsl/gsl_blas.h>
39
40namespace theplu {
41namespace yat {
42namespace utility {
43
44
45  Matrix::Matrix(void)
46    : blas_result_(NULL), m_(NULL)
47  {
48  }
49
50
51  Matrix::Matrix(const size_t& r, const size_t& c, double init_value)
52    : blas_result_(NULL), m_(NULL)
53  {
54    resize(r,c,init_value);
55  }
56
57
58  Matrix::Matrix(const Matrix& o)
59    : blas_result_(NULL), m_(o.create_gsl_matrix_copy())
60  {
61  }
62
63
64  // Constructor that gets data from istream
65  Matrix::Matrix(std::istream& is, char sep) 
66    throw (utility::IO_error,std::exception)
67    : blas_result_(NULL)
68  {
69    if (!is.good())
70      throw utility::IO_error("Matrix: istream is not good");
71
72    // read the data file and store in stl vectors (dynamically
73    // expandable)
74    std::vector<std::vector<double> > data_matrix;
75    unsigned int nof_columns=0;
76    unsigned 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      data_matrix.resize(data_matrix.size()+1);
85      std::vector<double>& v=data_matrix.back();
86      v.reserve(nof_columns);
87      std::string element;
88      std::stringstream ss(line);
89     
90      bool ok=true;
91      while(ok) {
92        if(sep=='\0')
93          ok=(ss>>element);
94        else
95          ok=getline(ss, element, sep);
96        if(!ok)
97          break;
98       
99        if (!element.size())
100          v.push_back(std::numeric_limits<double>::quiet_NaN());
101        else {
102          try {
103            v.push_back(convert<double>(element));
104          }
105          catch (std::runtime_error& e) {
106            std::stringstream ss(e.what());
107            ss << "\nMatrix.cc: " << element
108               << " is not accepted as a matrix element\n";
109            throw IO_error(ss.str());
110          }
111        }
112      }           
113      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
114          v.push_back(std::numeric_limits<double>::quiet_NaN());
115      if (!nof_columns)
116        nof_columns=v.size();
117      else if (v.size()!=nof_columns) {
118        std::ostringstream s;
119        s << "Matrix::Matrix(std::istream&, char) data file error: "
120          << "line " << nof_rows << " has " << v.size()
121          << " columns; expected " << nof_columns << " columns.";
122        throw utility::IO_error(s.str());
123      }
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(size_t i=0;i<nof_rows;i++)
141      for(size_t 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 column_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_column_iterator(&(*this)(0,i), this->columns());
183  }
184
185
186  Matrix::row_iterator Matrix::begin_row(size_t i)
187  {
188    return row_iterator(&(*this)(i,0), 1);
189  }
190
191
192  Matrix::const_row_iterator Matrix::begin_row(size_t i) const
193  {
194    return const_row_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    if (!m_)
220      return NULL;
221    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
222    if (!m)
223      throw utility::GSL_error("Matrix::create_gsl_matrix_copy failed to allocate memory");
224    if (gsl_matrix_memcpy(m,m_))
225      throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
226    return m;
227  }
228
229
230  void Matrix::delete_allocated_memory(void)
231  {
232    if (m_)
233      gsl_matrix_free(m_);
234    m_=NULL;
235  }
236
237
238  void Matrix::div(const Matrix& other)
239  {
240    assert(m_);
241    int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
242    if (status)
243      throw utility::GSL_error(std::string("matrix::div_elements",status));
244  }
245
246
247  Matrix::iterator Matrix::end(void)
248  {
249    return iterator(&(*this)(0,0)+rows()*columns(), 1);
250  }
251
252
253  Matrix::const_iterator Matrix::end(void) const
254  {
255    return const_iterator(&(*this)(0,0)+rows()*columns(), 1);
256  }
257
258
259  Matrix::column_iterator Matrix::end_column(size_t i)
260  {
261    return column_iterator(&(*this)(0,i)+rows()*columns(), this->columns());
262  }
263
264
265  Matrix::const_column_iterator Matrix::end_column(size_t i) const
266  {
267    return const_column_iterator(&(*this)(0,i)+rows()*columns(),this->columns());
268  }
269
270
271  Matrix::row_iterator Matrix::end_row(size_t i)
272  {
273    return row_iterator(&(*this)(i,0)+columns(), 1);
274  }
275
276
277  Matrix::const_row_iterator Matrix::end_row(size_t i) const
278  {
279    return const_row_iterator(&(*this)(i,0)+columns(), 1);
280  }
281
282
283  bool Matrix::equal(const Matrix& other, const double d) const
284  {
285    if (this==&other)
286      return true;
287    if (columns()!=other.columns() || rows()!=other.rows())
288      return false;
289    for (size_t i=0; i<rows(); i++)
290      if (!row_const_view(i).equal(other.row_const_view(i), d))
291        return false;
292    return true;
293  }
294
295
296  const gsl_matrix* Matrix::gsl_matrix_p(void) const
297  {
298    return m_;
299  }
300
301
302  gsl_matrix* Matrix::gsl_matrix_p(void)
303  {
304    return m_;
305  }
306
307
308  void Matrix::mul(const Matrix& other)
309  {
310    assert(m_);
311    int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());
312    if (status)
313      throw utility::GSL_error(std::string("Matrix::mul_elements",status));
314  }
315
316
317  void Matrix::resize(size_t r, size_t c, double init_value)
318  {
319    delete_allocated_memory();
320
321    assert( (r&&c) || (!r&&!c) );
322    if (r && c){
323      m_ = gsl_matrix_alloc(r,c);
324      if (!m_)
325        throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
326      all(init_value);
327    }
328
329    // no need to delete blas_result_ if the number of rows fit, it
330    // may be useful later.
331    if (blas_result_ && (blas_result_->size1!=rows())) {
332      gsl_matrix_free(blas_result_);
333      blas_result_=NULL;
334    }
335  }
336
337
338  size_t Matrix::rows(void) const
339  {
340    return (m_ ? m_->size1 : 0);
341  }
342
343
344  const VectorConstView Matrix::row_const_view(size_t col) const
345  {
346    return VectorConstView(*this, col, true);
347  }
348
349
350  VectorView Matrix::row_view(size_t row)
351  {
352    VectorView res(*this, row, true);
353    return res;
354  }
355
356
357  void Matrix::swap_columns(const size_t i, const size_t j)
358  {
359    assert(m_);
360    int status=gsl_matrix_swap_columns(m_, i, j);
361    if (status)
362      throw utility::GSL_error(std::string("Matrix::swap_columns",status));
363  }
364
365
366  void Matrix::swap_rowcol(const size_t i, const size_t j)
367  {
368    assert(m_);
369    int status=gsl_matrix_swap_rowcol(m_, i, j);
370    if (status)
371      throw utility::GSL_error(std::string("Matrix::swap_rowcol",status));
372  }
373
374
375  void Matrix::swap_rows(const size_t i, const size_t j)
376  {
377    assert(m_);
378    int status=gsl_matrix_swap_rows(m_, i, j);
379    if (status)
380      throw utility::GSL_error(std::string("Matrix::swap_rows",status));
381  }
382
383
384  void Matrix::transpose(void)
385  {
386    assert(m_);
387    if (columns()==rows())
388      gsl_matrix_transpose(m_); // this never fails
389    else {
390      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
391      if (!transposed)
392        throw utility::GSL_error("Matrix::transpose failed to allocate memory");
393      // next line never fails if allocation above succeeded.
394      gsl_matrix_transpose_memcpy(transposed,m_);
395      gsl_matrix_free(m_);
396      m_ = transposed;
397      if (blas_result_) {
398        gsl_matrix_free(blas_result_);
399        blas_result_=NULL;
400      }
401    }
402  }
403
404
405  double& Matrix::operator()(size_t row, size_t column)
406  {
407    assert(m_);
408    assert(row<rows());
409    assert(column<columns());
410    double* d=gsl_matrix_ptr(m_, row, column);
411    if (!d)
412      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
413    return *d;
414  }
415
416
417  const double& Matrix::operator()(size_t row, size_t column) const
418  {
419    assert(row<rows());
420    assert(column<columns());
421    const double* d=gsl_matrix_const_ptr(m_, row, column);
422    if (!d)
423      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
424    return *d;
425  }
426
427
428  bool Matrix::operator==(const Matrix& other) const
429  {
430    return equal(other);
431  }
432
433
434  bool Matrix::operator!=(const Matrix& other) const
435  {
436    return !equal(other);
437  }
438
439
440  const Matrix& Matrix::operator=( const Matrix& other )
441  {
442    if (this!=&other) {
443      if (rows()!=other.rows() || columns()!=other.columns())
444        resize(other.rows(), other.columns());
445       
446      if (m_)
447        if (gsl_matrix_memcpy(m_, other.gsl_matrix_p())) {
448          std::string s="Matrix::operator=  assignment failed";
449          throw utility::GSL_error(s);
450        }
451    }
452    return *this;
453  }
454
455
456  const Matrix& Matrix::operator+=(const Matrix& other)
457  {
458    assert(m_);
459    int status=gsl_matrix_add(m_, other.m_);
460    if (status)
461      throw utility::GSL_error(std::string("Matrix::operator+=", status));
462    return *this;
463  }
464
465
466  const Matrix& Matrix::operator+=(const double d)
467  {
468    assert(m_);
469    gsl_matrix_add_constant(m_, d);
470    return *this;
471  }
472
473
474  const Matrix& Matrix::operator-=(const Matrix& other)
475  {
476    assert(m_);
477    int status=gsl_matrix_sub(m_, other.m_);
478    if (status)
479      throw utility::GSL_error(std::string("Matrix::operator-=", status));
480    return *this;
481  }
482
483
484  const Matrix& Matrix::operator-=(const double d)
485  {
486    assert(m_);
487    gsl_matrix_add_constant(m_, -d);
488    return *this;
489  }
490
491
492  const Matrix& Matrix::operator*=(const Matrix& other)
493  {
494    assert(m_);
495    if ( blas_result_ && ((blas_result_->size1!=rows()) ||
496                          (blas_result_->size2!=other.columns())) ) {
497      gsl_matrix_free(blas_result_);
498      blas_result_=NULL;
499    }
500    if (!blas_result_) {
501      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
502      if (!blas_result_)
503        throw utility::GSL_error("Matrix::operator*= failed to allocate memory");
504    }
505    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, blas_result_);
506    gsl_matrix* tmp=m_;
507    m_ = blas_result_;
508    blas_result_=tmp;
509    return *this;
510  }
511
512
513  const Matrix& Matrix::operator*=(const double d)
514  {
515    assert(m_);
516    gsl_matrix_scale(m_, d);
517    return *this;
518  }
519
520
521  bool isnull(const Matrix& other)
522  {
523    return gsl_matrix_isnull(other.gsl_matrix_p());
524  }
525
526
527  double max(const Matrix& other)
528  {
529    return gsl_matrix_max(other.gsl_matrix_p());
530  }
531
532
533  double min(const Matrix& other)
534  {
535    return gsl_matrix_min(other.gsl_matrix_p());
536  }
537
538
539  void minmax_index(const Matrix& other,
540                    std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
541  {
542    gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,
543                            &max.first, &max.second);
544  }
545
546
547  bool nan(const Matrix& templat, Matrix& flag)
548  {
549    if (templat.rows()!=flag.rows() || templat.columns()!=flag.columns())
550      flag.resize(templat.rows(),templat.columns(),1.0);
551    return BinaryWeight()(templat.begin(), templat.end(), flag.begin());
552  }
553
554
555  void swap(Matrix& a, Matrix& b)
556  {
557    assert(a.gsl_matrix_p()); assert(b.gsl_matrix_p());
558    int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
559    if (status)
560      throw utility::GSL_error(std::string("swap(Matrix&,Matrix&)",status));
561  }
562
563
564  std::ostream& operator<<(std::ostream& s, const Matrix& m)
565  {
566    s.setf(std::ios::dec);
567    s.precision(12);
568    for(size_t i=0, j=0; i<m.rows(); i++)
569      for (j=0; j<m.columns(); j++) {
570        s << m(i,j);
571        if (j<m.columns()-1)
572          s << s.fill();
573        else if (i<m.rows()-1)
574          s << "\n";
575      }
576    return s;
577  }
578
579
580  Vector operator*(const Matrix& m, const VectorBase& v)
581  {
582    utility::Vector res(m.rows());
583    for (size_t i=0; i<res.size(); ++i)
584      res(i) = VectorConstView(m,i) * v;
585    return res;
586  }
587
588
589  Vector operator*(const VectorBase& v, const Matrix& m)
590  {
591    utility::Vector res(m.columns());
592    for (size_t i=0; i<res.size(); ++i)
593      res(i) = v * VectorConstView(m,i,false);
594    return res;
595  }
596
597}}} // of namespace utility, yat and thep
Note: See TracBrowser for help on using the repository browser.