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

Last change on this file since 2325 was 2325, checked in by Peter, 11 years ago

prefer using std::swap and avoid long lines

  • 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 2325 2010-09-22 15:05:37Z 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 <algorithm>
38#include <cassert>
39#include <climits>
40#include <cmath>
41#include <iostream>
42#include <limits>
43#include <sstream>
44#include <vector>
45
46namespace theplu {
47namespace yat {
48namespace utility {
49
50
51  Matrix::Matrix(void)
52    : blas_result_(NULL), m_(NULL)
53  {
54  }
55
56
57  Matrix::Matrix(const size_t& r, const size_t& c, double init_value)
58    : blas_result_(NULL), m_(NULL)
59  {
60    resize(r,c,init_value);
61  }
62
63
64  Matrix::Matrix(const Matrix& o)
65    : blas_result_(NULL), m_(o.create_gsl_matrix_copy())
66  {
67  }
68
69
70  // Constructor that gets data from istream
71  Matrix::Matrix(std::istream& is, char sep) 
72    throw (utility::IO_error,std::exception)
73    : blas_result_(NULL)
74  {
75    if (!is.good())
76      throw utility::IO_error("Matrix: istream is not good");
77
78    // read the data file and store in stl vectors (dynamically
79    // expandable)
80    std::vector<std::vector<double> > data_matrix;
81    try {
82      load(is, data_matrix, sep, '\n', true);
83    }
84    catch (utility::IO_error& e) {
85      std::stringstream ss(e.what());
86      ss << "\nMatrix(std::istream&): invalid dimensions\n";
87      throw IO_error(ss.str());
88    }
89    catch (runtime_error& e) {
90      std::stringstream ss(e.what());
91      ss << "\nMatrix(std::istream&): invalid matrix element\n";
92      throw IO_error(ss.str());
93    }
94
95    unsigned int nof_rows = data_matrix.size();
96    // if stream was empty, create nothing
97    if (!nof_rows)
98      return;
99
100    unsigned int nof_columns=data_matrix[0].size();
101
102    // convert the data to a gsl matrix
103    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
104    if (!m_)
105      throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
106
107    for (size_t i=0; i<data_matrix.size(); ++i) {
108      assert(data_matrix[i].size()==columns());
109      assert(i<rows());
110      std::copy(data_matrix[i].begin(), data_matrix[i].end(), begin_row(i));
111    }
112  }
113
114
115  Matrix::~Matrix(void)
116  {
117    delete_allocated_memory();
118    if (blas_result_)
119      gsl_matrix_free(blas_result_);
120    blas_result_=NULL;
121  }
122
123
124  void Matrix::all(const double value)
125  {
126    assert(m_);
127    gsl_matrix_set_all(m_, value);
128  }
129
130
131  Matrix::iterator Matrix::begin(void)
132  {
133    return iterator(&(*this)(0,0), 1);
134  }
135
136
137  Matrix::const_iterator Matrix::begin(void) const
138  {
139    return const_iterator(&(*this)(0,0), 1);
140  }
141
142
143  Matrix::column_iterator Matrix::begin_column(size_t i)
144  {
145    return column_iterator(&(*this)(0,i), this->columns());
146  }
147
148
149  Matrix::const_column_iterator Matrix::begin_column(size_t i) const
150  {
151    return const_column_iterator(&(*this)(0,i), this->columns());
152  }
153
154
155  Matrix::row_iterator Matrix::begin_row(size_t i)
156  {
157    return row_iterator(&(*this)(i,0), 1);
158  }
159
160
161  Matrix::const_row_iterator Matrix::begin_row(size_t i) const
162  {
163    return const_row_iterator(&(*this)(i,0), 1);
164  }
165
166
167  VectorView Matrix::column_view(size_t col)
168  {
169    VectorView res(*this, col, false);
170    return res;
171  }
172
173
174  const VectorConstView Matrix::column_const_view(size_t col) const
175  {
176    return VectorConstView(*this, col, false);
177  }
178
179
180  size_t Matrix::columns(void) const
181  {
182    return (m_ ? m_->size2 : 0);
183  }
184
185
186  gsl_matrix* Matrix::create_gsl_matrix_copy(void) const
187  {
188    if (!m_)
189      return NULL;
190    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
191    if (!m)
192      throw utility::GSL_error("Matrix::create_gsl_matrix_copy failed to allocate memory");
193    if (gsl_matrix_memcpy(m,m_))
194      throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
195    return m;
196  }
197
198
199  void Matrix::delete_allocated_memory(void)
200  {
201    if (m_)
202      gsl_matrix_free(m_);
203    m_=NULL;
204  }
205
206
207  void Matrix::div(const Matrix& other)
208  {
209    assert(m_);
210    int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
211    if (status)
212      throw utility::GSL_error(std::string("matrix::div_elements",status));
213  }
214
215
216  Matrix::iterator Matrix::end(void)
217  {
218    return iterator(&(*this)(0,0)+rows()*columns(), 1);
219  }
220
221
222  Matrix::const_iterator Matrix::end(void) const
223  {
224    return const_iterator(&(*this)(0,0)+rows()*columns(), 1);
225  }
226
227
228  Matrix::column_iterator Matrix::end_column(size_t i)
229  {
230    return column_iterator(&(*this)(0,i)+rows()*columns(), this->columns());
231  }
232
233
234  Matrix::const_column_iterator Matrix::end_column(size_t i) const
235  {
236    return const_column_iterator(&(*this)(0,i)+rows()*columns(),this->columns());
237  }
238
239
240  Matrix::row_iterator Matrix::end_row(size_t i)
241  {
242    return row_iterator(&(*this)(i,0)+columns(), 1);
243  }
244
245
246  Matrix::const_row_iterator Matrix::end_row(size_t i) const
247  {
248    return const_row_iterator(&(*this)(i,0)+columns(), 1);
249  }
250
251
252  bool Matrix::equal(const Matrix& other, const double d) const
253  {
254    if (this==&other)
255      return true;
256    if (columns()!=other.columns() || rows()!=other.rows())
257      return false;
258    for (size_t i=0; i<rows(); i++)
259      if (!row_const_view(i).equal(other.row_const_view(i), d))
260        return false;
261    return true;
262  }
263
264
265  const gsl_matrix* Matrix::gsl_matrix_p(void) const
266  {
267    return m_;
268  }
269
270
271  gsl_matrix* Matrix::gsl_matrix_p(void)
272  {
273    return m_;
274  }
275
276
277  void Matrix::mul(const Matrix& other)
278  {
279    assert(m_);
280    int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());
281    if (status)
282      throw utility::GSL_error(std::string("Matrix::mul_elements",status));
283  }
284
285
286  void Matrix::resize(size_t r, size_t c, double init_value)
287  {
288    delete_allocated_memory();
289
290    if (r && c){
291      m_ = gsl_matrix_alloc(r,c);
292      if (!m_)
293        throw utility::GSL_error("Matrix::resize failed to allocate memory");
294      all(init_value);
295    }
296    else if (r || c) {
297      std::stringstream msg;
298      msg << "Matrix::resize(" << r << ", " << c << "): incorrect dimensions\n";
299      throw runtime_error(msg.str());
300    }
301    // This should be triggered in else if above, right?
302    assert( (r&&c) || (!r&&!c) );
303
304    // no need to delete blas_result_ if the number of rows fit, it
305    // may be useful later.
306    if (blas_result_ && (blas_result_->size1!=rows())) {
307      gsl_matrix_free(blas_result_);
308      blas_result_=NULL;
309    }
310  }
311
312
313  size_t Matrix::rows(void) const
314  {
315    return (m_ ? m_->size1 : 0);
316  }
317
318
319  const VectorConstView Matrix::row_const_view(size_t col) const
320  {
321    return VectorConstView(*this, col, true);
322  }
323
324
325  VectorView Matrix::row_view(size_t row)
326  {
327    VectorView res(*this, row, true);
328    return res;
329  }
330
331
332  void Matrix::swap_columns(const size_t i, const size_t j)
333  {
334    assert(m_);
335    int status=gsl_matrix_swap_columns(m_, i, j);
336    if (status)
337      throw utility::GSL_error(std::string("Matrix::swap_columns",status));
338  }
339
340
341  void Matrix::swap_rowcol(const size_t i, const size_t j)
342  {
343    assert(m_);
344    int status=gsl_matrix_swap_rowcol(m_, i, j);
345    if (status)
346      throw utility::GSL_error(std::string("Matrix::swap_rowcol",status));
347  }
348
349
350  void Matrix::swap_rows(const size_t i, const size_t j)
351  {
352    assert(m_);
353    int status=gsl_matrix_swap_rows(m_, i, j);
354    if (status)
355      throw utility::GSL_error(std::string("Matrix::swap_rows",status));
356  }
357
358
359  void Matrix::transpose(void)
360  {
361    assert(m_);
362    if (columns()==rows())
363      gsl_matrix_transpose(m_); // this never fails
364    else {
365      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
366      if (!transposed)
367        throw utility::GSL_error("Matrix::transpose failed to allocate memory");
368      // next line never fails if allocation above succeeded.
369      gsl_matrix_transpose_memcpy(transposed,m_);
370      gsl_matrix_free(m_);
371      m_ = transposed;
372      if (blas_result_) {
373        gsl_matrix_free(blas_result_);
374        blas_result_=NULL;
375      }
376    }
377  }
378
379
380  double& Matrix::operator()(size_t row, size_t column)
381  {
382    assert(m_);
383    assert(row<rows());
384    assert(column<columns());
385    double* d=gsl_matrix_ptr(m_, row, column);
386    if (!d)
387      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
388    return *d;
389  }
390
391
392  const double& Matrix::operator()(size_t row, size_t column) const
393  {
394    assert(row<rows());
395    assert(column<columns());
396    const double* d=gsl_matrix_const_ptr(m_, row, column);
397    if (!d)
398      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
399    return *d;
400  }
401
402
403  bool Matrix::operator==(const Matrix& other) const
404  {
405    return equal(other);
406  }
407
408
409  bool Matrix::operator!=(const Matrix& other) const
410  {
411    return !equal(other);
412  }
413
414
415  const Matrix& Matrix::operator=( const Matrix& other )
416  {
417    if (this!=&other) {
418      if (rows()!=other.rows() || columns()!=other.columns())
419        resize(other.rows(), other.columns());
420       
421      if (m_)
422        if (gsl_matrix_memcpy(m_, other.gsl_matrix_p())) {
423          std::string s="Matrix::operator=  assignment failed";
424          throw utility::GSL_error(s);
425        }
426    }
427    return *this;
428  }
429
430
431  const Matrix& Matrix::operator+=(const Matrix& other)
432  {
433    assert(m_);
434    int status=gsl_matrix_add(m_, other.m_);
435    if (status)
436      throw utility::GSL_error(std::string("Matrix::operator+=", status));
437    return *this;
438  }
439
440
441  const Matrix& Matrix::operator+=(const double d)
442  {
443    assert(m_);
444    gsl_matrix_add_constant(m_, d);
445    return *this;
446  }
447
448
449  const Matrix& Matrix::operator-=(const Matrix& other)
450  {
451    assert(m_);
452    int status=gsl_matrix_sub(m_, other.m_);
453    if (status)
454      throw utility::GSL_error(std::string("Matrix::operator-=", status));
455    return *this;
456  }
457
458
459  const Matrix& Matrix::operator-=(const double d)
460  {
461    assert(m_);
462    gsl_matrix_add_constant(m_, -d);
463    return *this;
464  }
465
466
467  const Matrix& Matrix::operator*=(const Matrix& other)
468  {
469    assert(m_);
470    if ( blas_result_ && ((blas_result_->size1!=rows()) ||
471                          (blas_result_->size2!=other.columns())) ) {
472      gsl_matrix_free(blas_result_);
473      blas_result_=NULL;
474    }
475    if (!blas_result_) {
476      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
477      if (!blas_result_)
478        throw utility::GSL_error("Matrix::operator*= failed to allocate memory");
479    }
480    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, 
481                   blas_result_);
482    std::swap(m_, blas_result_);
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.