My Project
gsl_matrix.hh
Go to the documentation of this file.
1 /* -*- mia-c++ -*-
2  *
3  * This file is part of MIA - a toolbox for medical image analysis
4  * Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
5  *
6  * MIA is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with MIA; if not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef GSLPP_MATRIX_HH
22 #define GSLPP_MATRIX_HH
23 
24 #include <cassert>
25 #include <iterator>
26 #include <gsl/gsl_matrix.h>
27 #include <mia/core/gsl_vector.hh>
28 
29 namespace gsl
30 {
31 
33 {
34  friend class const_matrix_iterator;
35 public:
36  matrix_iterator(gsl_matrix *m, bool begin):
37  m_matrix(m),
38  m_current(begin ? m->data : m->data + m->size1 * m->tda),
39  m_current_column(0),
40  m_row_jump(m->tda - m->size2)
41  {
42  }
43 
45  m_matrix(nullptr),
46  m_current(nullptr),
47  m_current_column(0)
48  {
49  }
50 
51 
52  matrix_iterator(const matrix_iterator& other) = default;
53 
54  double& operator *()
55  {
56  assert(m_current);
57  return *m_current;
58  }
59 
60  matrix_iterator& operator ++()
61  {
62  assert(m_current);
63  ++m_current;
64  ++m_current_column;
65 
66  if (m_current_column == m_matrix->size2) {
67  m_current += m_row_jump;
68  m_current_column = 0;
69  }
70 
71  return *this;
72  }
73 
74 
75  matrix_iterator operator ++(int)
76  {
77  matrix_iterator result(*this);
78  ++(*this);
79  return result;
80  }
81 
82  friend bool operator == (const matrix_iterator& lhs, const matrix_iterator& rhs);
83  friend bool operator != (const matrix_iterator& lhs, const matrix_iterator& rhs);
84 private:
85  gsl_matrix *m_matrix;
86  double *m_current;
87  size_t m_current_column;
88  size_t m_row_jump;
89 };
90 
93 
94 
96 {
97 public:
98  const_matrix_iterator(const gsl_matrix *m, bool begin):
99  m_matrix(m),
100  m_current(begin ? m->data : m->data + m->size1 * m->tda),
101  m_current_column(0),
102  m_row_jump(m->tda - m->size2)
103  {
104  }
105 
107  m_matrix(other.m_matrix),
108  m_current(other.m_current),
109  m_current_column(other.m_current_column),
110  m_row_jump(other.m_row_jump)
111  {
112  }
113 
115  m_matrix(nullptr),
116  m_current(nullptr),
117  m_current_column(0)
118  {
119  }
120 
121 
123 
124  double operator *() const
125  {
126  assert(m_current);
127  return *m_current;
128  }
129 
130  const_matrix_iterator& operator ++()
131  {
132  assert(m_current);
133  ++m_current;
134  ++m_current_column;
135 
136  if (m_current_column == m_matrix->size2) {
137  m_current += m_row_jump;
138  m_current_column = 0;
139  }
140 
141  return *this;
142  }
143 
144  const_matrix_iterator operator ++(int)
145  {
146  const_matrix_iterator result(*this);
147  ++(*this);
148  return result;
149  }
150 
151  friend bool operator == (const const_matrix_iterator& lhs, const const_matrix_iterator& rhs);
152  friend bool operator != (const const_matrix_iterator& lhs, const const_matrix_iterator& rhs);
153 private:
154  const gsl_matrix *m_matrix;
155  const double *m_current;
156  size_t m_current_column;
157  size_t m_row_jump;
158 };
159 
162 
163 
164 
165 
166 
172 {
173 public:
176 
178 
185  Matrix(size_t rows, size_t columns, bool clean);
186 
193  Matrix(size_t rows, size_t columns, double init);
194 
202  Matrix(size_t rows, size_t columns, const double *init);
203 
207  Matrix(const Matrix& other);
208 
209 
213  Matrix(Matrix&& other);
214 
215 
220  Matrix(gsl_matrix *m);
221 
226  Matrix(const gsl_matrix *m);
227 
233 
234 
238  Matrix& operator =(const Matrix& other);
239 
243  Matrix& operator =(Matrix&& other);
244 
245 
252  void reset(size_t rows, size_t columns, bool clean);
253 
260  void reset(size_t rows, size_t columns, double init);
261 
263 
264 
268  size_t rows()const;
269 
273  size_t cols()const;
274 
281  void set(size_t i, size_t j, double x);
282 
290  double operator ()(size_t i, size_t j) const;
291 
295  operator gsl_matrix *();
296 
300  operator const gsl_matrix *() const;
301 
307 
313 
320 
327 
334 
341 
348  void set_row(int r, const Vector& row);
349 
356 
362  ConstVectorView get_row(int r) const;
363 
370  void set_column(int c, const Vector& col);
371 
378 
385 
393  double dot_row(int r, const Vector& v) const;
394 
403  double dot_column(int c, const Vector& col) const;
404 
409  void print(std::ostream& os) const;
410 
416  Matrix& operator -=(const Matrix& rhs)
417  {
418  gsl_matrix_sub(*this, rhs);
419  return *this;
420  }
421 
427  Matrix& operator +=(const Matrix& rhs)
428  {
429  gsl_matrix_add(*this, rhs);
430  return *this;
431  }
432 
438  Matrix& operator *=(double rhs)
439  {
440  gsl_matrix_scale(*this, rhs);
441  return *this;
442  }
443 
444  bool is_valid() const;
445 
446  bool is_writable() const;
447 
448 private:
449  gsl_matrix *m_matrix;
450  const gsl_matrix *m_const_matrix;
451  bool m_owner;
452 };
453 
454 
455 Matrix EXPORT_GSL operator * (const Matrix& lhs, const Matrix& rhs);
456 Matrix EXPORT_GSL operator + (const Matrix& lhs, const Matrix& rhs);
457 Matrix EXPORT_GSL operator - (const Matrix& lhs, const Matrix& rhs);
458 
459 inline std::ostream& operator << (std::ostream& os, const Matrix& m)
460 {
461  m.print(os);
462  return os;
463 }
464 
470 
473 };
474 
480 
481 
482 } // end namespace
483 
484 namespace std
485 {
486 
487 template <>
488 class iterator_traits< gsl::const_matrix_iterator >
489 {
490 public:
491  typedef size_t difference_type;
492  typedef double value_type;
493  typedef const double *pointer;
494  typedef const double& reference;
495  typedef forward_iterator_tag iterator_category;
496 };
497 
498 template <>
499 class iterator_traits< gsl::matrix_iterator >
500 {
501 public:
502  typedef size_t difference_type;
503  typedef double value_type;
504  typedef double *pointer;
505  typedef double& reference;
506  typedef forward_iterator_tag iterator_category;
507 };
508 
509 }
510 
511 
512 #endif
EXPORT_2D C2DFVectorfield & operator+=(C2DFVectorfield &a, const C2DFVectorfield &b)
size_t rows() const
const_matrix_iterator begin() const
Matrix(gsl_matrix *m)
double dot_row(int r, const Vector &v) const
Matrix row_covariance() const
Matrix transposed() const
const_matrix_iterator const_iterator
Definition: gsl_matrix.hh:175
VectorView get_column(int c)
void set_row(int r, const Vector &row)
void reset(size_t rows, size_t columns, double init)
ConstVectorView get_column(int c) const
VectorView get_row(int r)
Matrix column_covariance() const
Matrix(size_t rows, size_t columns, const double *init)
const_matrix_iterator end() const
matrix_iterator iterator
Definition: gsl_matrix.hh:174
void set_column(int c, const Vector &col)
Matrix(Matrix &&other)
Matrix(const gsl_matrix *m)
void print(std::ostream &os) const
double dot_column(int c, const Vector &col) const
matrix_iterator end()
Matrix(const Matrix &other)
ConstVectorView get_row(int r) const
void reset(size_t rows, size_t columns, bool clean)
matrix_iterator begin()
Matrix(size_t rows, size_t columns, bool clean)
bool is_valid() const
void set(size_t i, size_t j, double x)
size_t cols() const
bool is_writable() const
Matrix(size_t rows, size_t columns, double init)
const_matrix_iterator(const gsl_matrix *m, bool begin)
Definition: gsl_matrix.hh:98
const_matrix_iterator(const const_matrix_iterator &other)=default
const_matrix_iterator(const matrix_iterator &other)
Definition: gsl_matrix.hh:106
matrix_iterator(const matrix_iterator &other)=default
matrix_iterator(gsl_matrix *m, bool begin)
Definition: gsl_matrix.hh:36
#define EXPORT_GSL
Definition: gsl_defines.hh:38
std::ostream & operator<<(std::ostream &os, const Matrix &m)
Definition: gsl_matrix.hh:459
bool EXPORT_GSL operator==(const matrix_iterator &lhs, const matrix_iterator &rhs)
bool EXPORT_GSL operator!=(const matrix_iterator &lhs, const matrix_iterator &rhs)
vector_iterator operator-(const vector_iterator &it, int dist)
vector_iterator operator+(const vector_iterator &it, int dist)
Matrix EXPORT_GSL operator*(const Matrix &lhs, const Matrix &rhs)
void EXPORT_GSL matrix_inv_sqrt(Matrix &m)
CSymmvEvalEvec(Matrix m)