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>
28
29namespace gsl
30{
31
33{
35public:
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);
84private:
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{
97public:
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
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);
153private:
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{
173public:
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
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
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
448private:
449 gsl_matrix *m_matrix;
450 const gsl_matrix *m_const_matrix;
451 bool m_owner;
452};
453
454
455Matrix EXPORT_GSL operator * (const Matrix& lhs, const Matrix& rhs);
456Matrix EXPORT_GSL operator + (const Matrix& lhs, const Matrix& rhs);
457Matrix EXPORT_GSL operator - (const Matrix& lhs, const Matrix& rhs);
458
459inline std::ostream& operator << (std::ostream& os, const Matrix& m)
460{
461 m.print(os);
462 return os;
463}
464
474
480
481
482} // end namespace
483
484namespace std
485{
486
487template <>
488class iterator_traits< gsl::const_matrix_iterator >
489{
490public:
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
498template <>
499class iterator_traits< gsl::matrix_iterator >
500{
501public:
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
bool operator!=(const C2DImage &a, const C2DImage &b)
Definition 2d/image.hh:413
EXPORT_2D C2DFVectorfield & operator+=(C2DFVectorfield &a, const C2DFVectorfield &b)
bool operator==(const CAttribute &a, const CAttribute &b)
Definition attributes.hh:95
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
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
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)
matrix_iterator(const matrix_iterator &other)=default
matrix_iterator(gsl_matrix *m, bool begin)
Definition gsl_matrix.hh:36
#define EXPORT_GSL
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)
std::ostream & operator<<(std::ostream &os, const Matrix &m)
void EXPORT_GSL matrix_inv_sqrt(Matrix &m)
F operator*(const typename TSparseSolver< F >::A_mult_x &A, const F &x)
CSymmvEvalEvec(Matrix m)