1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
|
#ifndef _NL_MATRIX_H_
#define _NL_MATRIX_H_
#include <cassert>
#include <utility>
#include <algorithm>
#include <gsl/gsl_matrix.h>
namespace Geom { namespace NL {
class Matrix;
void swap(Matrix & m1, Matrix & m2);
class Matrix
{
public:
// the matrix is not inizialized
Matrix(size_t n1, size_t n2)
: m_rows(n1), m_columns(n2)
{
m_matrix = gsl_matrix_alloc(n1, n2);
}
Matrix(size_t n1, size_t n2, double x)
:m_rows(n1), m_columns(n2)
{
m_matrix = gsl_matrix_alloc(n1, n2);
gsl_matrix_set_all(m_matrix, x);
}
Matrix( Matrix const& _matrix )
: m_rows(_matrix.rows()), m_columns(_matrix.columns())
{
m_matrix = gsl_matrix_alloc(rows(), columns());
gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
}
// explicit
// Matrix( gsl_matrix* m, size_t n1, size_t n2)
// : m_rows(n1), m_columns(n2), m_matrix(m)
// {
// }
Matrix & operator=(Matrix const& _matrix)
{
assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
return *this;
}
virtual ~Matrix()
{
gsl_matrix_free(m_matrix);
}
void set_all( double x = 0)
{
gsl_matrix_set_all(m_matrix, x);
}
void set_identity()
{
gsl_matrix_set_identity(m_matrix);
}
const double & operator() (size_t i, size_t j) const
{
return *gsl_matrix_const_ptr(m_matrix, i, j);
}
double & operator() (size_t i, size_t j)
{
return *gsl_matrix_ptr(m_matrix, i, j);
}
gsl_matrix* get_gsl_matrix()
{
return m_matrix;
}
void swap_rows(size_t i, size_t j)
{
gsl_matrix_swap_rows(m_matrix, i, j);
}
void swap_columns(size_t i, size_t j)
{
gsl_matrix_swap_columns(m_matrix, i, j);
}
Matrix & transpose()
{
gsl_matrix_transpose(m_matrix);
return (*this);
}
Matrix & scale(double x)
{
gsl_matrix_scale(m_matrix, x);
return (*this);
}
Matrix & translate(double x)
{
gsl_matrix_add_constant(m_matrix, x);
return (*this);
}
Matrix & operator+=(Matrix const& _matrix)
{
gsl_matrix_add(m_matrix, _matrix.m_matrix);
return (*this);
}
Matrix & operator-=(Matrix const& _matrix)
{
gsl_matrix_sub(m_matrix, _matrix.m_matrix);
return (*this);
}
bool is_zero() const
{
return gsl_matrix_isnull(m_matrix);
}
bool is_positive() const
{
return gsl_matrix_ispos(m_matrix);
}
bool is_negative() const
{
return gsl_matrix_isneg(m_matrix);
}
bool is_non_negative() const
{
for ( unsigned int i = 0; i < rows(); ++i )
{
for ( unsigned int j = 0; j < columns(); ++j )
{
if ( (*this)(i,j) < 0 ) return false;
}
}
return true;
}
double max() const
{
return gsl_matrix_max(m_matrix);
}
double min() const
{
return gsl_matrix_min(m_matrix);
}
std::pair<size_t, size_t>
max_index() const
{
std::pair<size_t, size_t> indices;
gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));
return indices;
}
std::pair<size_t, size_t>
min_index() const
{
std::pair<size_t, size_t> indices;
gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));
return indices;
}
friend
void swap(Matrix & m1, Matrix & m2);
size_t rows() const
{
return m_rows;
}
size_t columns() const
{
return m_columns;
}
private:
size_t m_rows, m_columns;
gsl_matrix* m_matrix;
};
void swap(Matrix & m1, Matrix & m2)
{
assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );
std::swap(m1.m_matrix, m2.m_matrix);
}
} } // end namespaces
#endif /*_NL_MATRIX_H_*/
|
