source: trunk/yat/utility/Matrix.cc

Last change on this file was 4060, checked in by Peter, 6 weeks ago

prefer throw_with_nested, when info in the caught expression makes sense to expose

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