yat  0.16pre
Matrix.h
1 #ifndef _theplu_yat_utility_matrix_
2 #define _theplu_yat_utility_matrix_
3 
4 // $Id: Matrix.h 3735 2018-05-08 01:45:58Z peter $
5 
6 /*
7  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
8  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
10  Copyright (C) 2007, 2008, 2009 Jari Häkkinen, Peter Johansson
11  Copyright (C) 2012, 2017 Peter Johansson
12 
13  This file is part of the yat library, http://dev.thep.lu.se/yat
14 
15  The yat library is free software; you can redistribute it and/or
16  modify it under the terms of the GNU General Public License as
17  published by the Free Software Foundation; either version 3 of the
18  License, or (at your option) any later version.
19 
20  The yat library is distributed in the hope that it will be useful,
21  but WITHOUT ANY WARRANTY; without even the implied warranty of
22  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
23  General Public License for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with yat. If not, see <http://www.gnu.org/licenses/>.
27 */
28 
29 #include <cassert>
30 
31 #include "config_public.h"
32 
33 // include these operations as they are expected when using Matrix
34 #include "BLAS_level2.h"
35 #include "BLAS_level3.h"
36 
37 #include "Exception.h"
38 #include "MatrixExpression.h"
39 #include "StrideIterator.h"
40 #include "Vector.h"
41 #include "VectorConstView.h"
42 #include "VectorExpression.h"
43 #include "VectorView.h"
44 #include "yat_assert.h"
45 
46 #include <gsl/gsl_matrix.h>
47 
48 #include <cstddef> // size_t
49 #include <iosfwd>
50 
51 namespace theplu {
52 namespace yat {
53 namespace utility {
54 
55  class VectorBase;
56 
74  class Matrix : public BasicMatrix<Matrix>
75  {
76  public:
82  typedef double value_type;
83 
89  typedef double& reference;
90 
96  typedef const double& const_reference;
97 
102 
107 
112 
117 
122 
127 
134  Matrix(void);
135 
144  Matrix(const size_t& r, const size_t& c, double init_value=0);
145 
152  Matrix(const Matrix&);
153 
172  template<class T>
173  Matrix(const MatrixExpression<T>& other);
174 
175 #ifdef YAT_HAVE_RVALUE
176 
181  Matrix(Matrix&&) noexcept;
182 
198  template<class T>
199  Matrix(MatrixExpression<T>&& other)
200  : blas_result_(NULL), m_(NULL)
201  {
202  other.move(m_);
203  }
204 
205 #endif
206 
237  explicit Matrix(std::istream &, char sep='\0')
238  throw(utility::IO_error, std::exception);
239 
243  ~Matrix(void);
244 
248  void all(const double value);
249 
256  iterator begin(void);
257 
264  const_iterator begin(void) const;
265 
271  iterator begin_column(size_t i);
272 
278  const_iterator begin_column(size_t i) const;
279 
285  iterator begin_row(size_t i);
286 
292  const_iterator begin_row(size_t i) const;
293 
297  VectorView column_view(size_t i);
298 
302  const VectorConstView column_const_view(size_t) const;
303 
307  size_t columns(void) const;
308 
316  void div(const Matrix& b);
317 
321  iterator end(void);
322 
326  const_iterator end(void) const;
327 
331  iterator end_column(size_t i);
332 
336  const_iterator end_column(size_t i) const;
337 
341  iterator end_row(size_t i);
342 
346  const_iterator end_row(size_t i) const;
347 
357  bool equal(const Matrix&, const double precision=0) const;
358 
362  const gsl_matrix* gsl_matrix_p(void) const;
363 
367  gsl_matrix* gsl_matrix_p(void);
368 
376  void mul(const Matrix& b);
377 
388  void resize(size_t r, size_t c, double init_value=0);
389 
393  size_t rows(void) const;
394 
398  VectorView row_view(size_t);
399 
403  const VectorConstView row_const_view(size_t) const;
404 
410  void swap_columns(const size_t i,const size_t j);
411 
420  void swap_rowcol(const size_t i,const size_t j);
421 
427  void swap_rows(const size_t i, const size_t j);
428 
438  void transpose(void);
439 
449  double& operator()(size_t row,size_t column);
450 
461  const double& operator()(size_t row,size_t column) const;
462 
474  bool operator==(const Matrix& other) const;
475 
487  bool operator!=(const Matrix& other) const;
488 
497  const Matrix& operator=(const Matrix& other);
498 
517  template<class T>
518  Matrix& operator=(const MatrixExpression<T>& rhs);
519 
520 #ifdef YAT_HAVE_RVALUE
521 
526  Matrix& operator=(Matrix&& other);
527 
546  template<class T>
547  Matrix& operator=(MatrixExpression<T>&& rhs)
548  {
549  // This function steals data from rhs into this, and give the old
550  // m_ in return to rhs so destructor of rhs takes care of
551  // deallocation.
552 
553  rhs.move(m_);
554  return *this;
555  }
556 
557 #endif
558 
570  const Matrix& operator+=(const Matrix& b);
571 
577  template<class T>
578  Matrix& operator+=(const MatrixExpression<T>& rhs)
579  {
580  *this = *this + rhs;
581  return *this;
582  }
583 
590  const Matrix& operator+=(const double d);
591 
603  const Matrix& operator-=(const Matrix&);
604 
610  template<class T>
611  Matrix& operator-=(const MatrixExpression<T>& rhs)
612  {
613  *this = *this - rhs;
614  return *this;
615  }
616 
623  const Matrix& operator-=(const double d);
624 
634  const Matrix& operator*=(const Matrix&);
635 
644  const Matrix& operator*=(double d);
645 
646  private:
647 
659  gsl_matrix* create_gsl_matrix_copy(void) const;
660 
666  void delete_allocated_memory(void);
667 
668  // blas_result_ is used to temporarily store result in BLAS calls
669  // and when not NULL it should always have the same size as
670  // m_. Memory is not allocated for blas_result_ until it is
671  // needed.
672  gsl_matrix* blas_result_;
673  gsl_matrix* m_;
674  };
675 
683  void inverse_svd(const Matrix& input, Matrix& result);
684 
693  bool isnull(const Matrix&);
694 
702  double max(const Matrix&);
703 
711  double min(const Matrix&);
712 
723  void minmax_index(const Matrix&,
724  std::pair<size_t,size_t>& min,
725  std::pair<size_t,size_t>& max);
726 
742  bool nan(const Matrix& templat, Matrix& flag);
743 
757  void swap(Matrix&, Matrix&);
758 
764  std::ostream& operator<< (std::ostream& s, const Matrix&);
765 
775  double trace(const Matrix&);
776 
778  template<class T>
780  : blas_result_(NULL), m_(NULL)
781  {
782  // In principle, we could implement this as a move constructor
783  // and steal the data (rather than copy), but it's a bit hackish
784  // to workaround the constness, so users who are eager for speed
785  // should enable rvalues and get access to the faster move
786  // constructor below.
787 
788  try {
789  detail::copy(m_, other.gsl_matrix_p());
790  }
791  catch (GSL_error& e) {
792  detail::deallocate(m_);
793  throw e;
794  }
795  }
796 
797 
798  template<class T>
800  {
801  // Be careful because *this might be embedded in rhs.
802  rhs.get(blas_result_);
803  std::swap(m_, blas_result_);
804  return *this;
805  }
806 
807 
808 }}} // of namespace utility, yat, and theplu
809 
810 #endif
const VectorConstView column_const_view(size_t) const
double trace(const Matrix &)
Trace of matrix.
void mul(const Matrix &b)
const double & const_reference
Definition: Matrix.h:96
void inverse_svd(const Matrix &input, Matrix &result)
iterator end_row(size_t i)
Class for errors reported from underlying GSL calls.
Definition: Exception.h:87
The Department of Theoretical Physics namespace as we define it.
StrideIterator< const double * > const_iterator
Definition: Matrix.h:106
This is the yat interface to gsl_vector_view.
Definition: VectorView.h:79
Definition: stl_utility.h:64
void resize(size_t r, size_t c, double init_value=0)
Resize Matrix.
size_t rows(void) const
iterator begin_row(size_t i)
const gsl_matrix * gsl_matrix_p(void) const
Definition: MatrixExpression.h:140
StrideIterator< const double * > const_row_iterator
Definition: Matrix.h:126
bool nan(const Matrix &templat, Matrix &flag)
Create a Matrix flag indicating NaN&#39;s in another Matrix templat.
StrideIterator< double * > row_iterator
Definition: Matrix.h:121
double & reference
Definition: Matrix.h:89
void minmax_index(const Matrix &, std::pair< size_t, size_t > &min, std::pair< size_t, size_t > &max)
Locate the maximum and minimum element in the Matrix.
Matrix(void)
The default constructor.
void all(const double value)
Matrix & operator-=(const MatrixExpression< T > &rhs)
Subtraction and assign operator.
Definition: Matrix.h:611
Matrix & operator+=(const MatrixExpression< T > &rhs)
Addition and assign operator.
Definition: Matrix.h:578
Read-only view.
Definition: VectorConstView.h:56
iterator begin_column(size_t i)
double min(const Matrix &)
Get the minimum value of the Matrix.
double value_type
Definition: Matrix.h:82
bool equal(const Matrix &, const double precision=0) const
Check whether matrices are equal within a user defined precision, set by precision.
void swap(Matrix &, Matrix &)
Exchange all elements between the matrices by copying.
StrideIterator< const double * > const_column_iterator
Definition: Matrix.h:116
const Matrix & operator-=(const Matrix &)
Subtract and assign operator.
std::ostream & operator<<(std::ostream &s, const Matrix &)
The output operator for the Matrix class.
VectorView column_view(size_t i)
An expression that can be converted to a Matrix.
Definition: MatrixExpression.h:46
void get(gsl_matrix *&m) const
Definition: MatrixExpression.h:112
void div(const Matrix &b)
Class to report errors associated with IO operations.
Definition: Exception.h:109
VectorView row_view(size_t)
const VectorConstView row_const_view(size_t) const
Interface to GSL matrix.
Definition: Matrix.h:74
bool isnull(const Matrix &)
Check if all elements of the Matrix are zero.
Definition: BasicMatrix.h:38
const Matrix & operator=(const Matrix &other)
The assignment operator.
const gsl_matrix * gsl_matrix_p(void) const
StrideIterator< double * > iterator
Definition: Matrix.h:101
StrideIterator< double * > column_iterator
Definition: Matrix.h:111
const Matrix & operator+=(const Matrix &b)
Add and assign operator.
void swap_rows(const size_t i, const size_t j)
Swap rows i and j.
const Matrix & operator*=(const Matrix &)
Multiply and assignment operator.
iterator end_column(size_t i)
void swap_columns(const size_t i, const size_t j)
Swap columns i and j.
void transpose(void)
Transpose the matrix.
Adaptor using a stride on underlying iterator.
Definition: StrideIterator.h:50
double max(const Matrix &)
Get the maximum value of the Matrix.
void swap_rowcol(const size_t i, const size_t j)
Swap row i and column j.
size_t columns(void) const

Generated on Thu Jan 10 2019 02:54:19 for yat by  doxygen 1.8.11