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

Last change on this file since 2748 was 2748, checked in by Peter, 9 years ago

avoid extra if statement

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