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

Last change on this file since 1700 was 1700, checked in by Peter, 13 years ago

reserve size in vector to avoid copying

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