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

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

fixes #281. Change all throws of std::runtime_error to theplu::yat::utility::runtime_error to clarify that the error comes from yat. Also removed some throw declarations.

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