1 | /*
|
---|
2 | * MatrixContent.cpp
|
---|
3 | *
|
---|
4 | * Created on: Nov 14, 2010
|
---|
5 | * Author: heber
|
---|
6 | */
|
---|
7 |
|
---|
8 |
|
---|
9 | // include config.h
|
---|
10 | #ifdef HAVE_CONFIG_H
|
---|
11 | #include <config.h>
|
---|
12 | #endif
|
---|
13 |
|
---|
14 | #include "Helpers/MemDebug.hpp"
|
---|
15 |
|
---|
16 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
17 | #include "Helpers/Assert.hpp"
|
---|
18 | #include "Exceptions/NotInvertibleException.hpp"
|
---|
19 | #include "Helpers/fast_functions.hpp"
|
---|
20 | #include "Helpers/Assert.hpp"
|
---|
21 | #include "LinearAlgebra/Vector.hpp"
|
---|
22 | #include "LinearAlgebra/VectorContent.hpp"
|
---|
23 | #include "LinearAlgebra/MatrixContent.hpp"
|
---|
24 |
|
---|
25 | #include <gsl/gsl_blas.h>
|
---|
26 | #include <gsl/gsl_eigen.h>
|
---|
27 | #include <gsl/gsl_matrix.h>
|
---|
28 | #include <gsl/gsl_multimin.h>
|
---|
29 | #include <gsl/gsl_vector.h>
|
---|
30 | #include <cmath>
|
---|
31 | #include <iostream>
|
---|
32 |
|
---|
33 | using namespace std;
|
---|
34 |
|
---|
35 |
|
---|
36 | /** Constructor for class MatrixContent.
|
---|
37 | * \param rows number of rows
|
---|
38 | * \param columns number of columns
|
---|
39 | */
|
---|
40 | MatrixContent::MatrixContent(size_t _rows, size_t _columns) :
|
---|
41 | rows(_rows),
|
---|
42 | columns(_columns)
|
---|
43 | {
|
---|
44 | content = gsl_matrix_calloc(rows, columns);
|
---|
45 | }
|
---|
46 |
|
---|
47 | /** Constructor for class MatrixContent.
|
---|
48 | * \param rows number of rows
|
---|
49 | * \param columns number of columns
|
---|
50 | * \param *src array with components to initialize matrix with
|
---|
51 | */
|
---|
52 | MatrixContent::MatrixContent(size_t _rows, size_t _columns, const double *src) :
|
---|
53 | rows(_rows),
|
---|
54 | columns(_columns)
|
---|
55 | {
|
---|
56 | content = gsl_matrix_calloc(rows, columns);
|
---|
57 | set(0,0, src[0]);
|
---|
58 | set(1,0, src[1]);
|
---|
59 | set(2,0, src[2]);
|
---|
60 |
|
---|
61 | set(0,1, src[3]);
|
---|
62 | set(1,1, src[4]);
|
---|
63 | set(2,1, src[5]);
|
---|
64 |
|
---|
65 | set(0,2, src[6]);
|
---|
66 | set(1,2, src[7]);
|
---|
67 | set(2,2, src[8]);
|
---|
68 | }
|
---|
69 |
|
---|
70 | /** Constructor for class MatrixContent.
|
---|
71 | * We embed the given gls_matrix pointer within this class and set it to NULL
|
---|
72 | * afterwards.
|
---|
73 | * \param *src source gsl_matrix vector to embed within this class
|
---|
74 | */
|
---|
75 | MatrixContent::MatrixContent(gsl_matrix *&src) :
|
---|
76 | rows(src->size1),
|
---|
77 | columns(src->size2)
|
---|
78 | {
|
---|
79 | content = gsl_matrix_alloc(src->size1, src->size2);
|
---|
80 | gsl_matrix_memcpy(content,src);
|
---|
81 | // content = src;
|
---|
82 | // src = NULL;
|
---|
83 | }
|
---|
84 |
|
---|
85 | /** Copy constructor for class MatrixContent.
|
---|
86 | * \param &src reference to source MatrixContent
|
---|
87 | */
|
---|
88 | MatrixContent::MatrixContent(const MatrixContent &src) :
|
---|
89 | rows(src.rows),
|
---|
90 | columns(src.columns)
|
---|
91 | {
|
---|
92 | content = gsl_matrix_alloc(src.rows, src.columns);
|
---|
93 | gsl_matrix_memcpy(content,src.content);
|
---|
94 | }
|
---|
95 |
|
---|
96 | /** Copy constructor for class MatrixContent.
|
---|
97 | * \param *src pointer to source MatrixContent
|
---|
98 | */
|
---|
99 | MatrixContent::MatrixContent(const MatrixContent *src) :
|
---|
100 | rows(src->rows),
|
---|
101 | columns(src->columns)
|
---|
102 | {
|
---|
103 | ASSERT(src != NULL, "MatrixContent::MatrixContent - pointer to source matrix is NULL!");
|
---|
104 | content = gsl_matrix_alloc(src->rows, src->columns);
|
---|
105 | gsl_matrix_memcpy(content,src->content);
|
---|
106 | }
|
---|
107 |
|
---|
108 | /** Destructor for class MatrixContent.
|
---|
109 | */
|
---|
110 | MatrixContent::~MatrixContent()
|
---|
111 | {
|
---|
112 | gsl_matrix_free(content);
|
---|
113 | }
|
---|
114 |
|
---|
115 | /** Set matrix to identity.
|
---|
116 | */
|
---|
117 | void MatrixContent::setIdentity()
|
---|
118 | {
|
---|
119 | for(int i=rows;i--;){
|
---|
120 | for(int j=columns;j--;){
|
---|
121 | set(i,j,i==j);
|
---|
122 | }
|
---|
123 | }
|
---|
124 | }
|
---|
125 |
|
---|
126 | /** Set all matrix components to zero.
|
---|
127 | */
|
---|
128 | void MatrixContent::setZero()
|
---|
129 | {
|
---|
130 | for(int i=rows;i--;){
|
---|
131 | for(int j=columns;j--;){
|
---|
132 | set(i,j,0.);
|
---|
133 | }
|
---|
134 | }
|
---|
135 | }
|
---|
136 |
|
---|
137 | /** Copy operator for MatrixContent with self-assignment check.
|
---|
138 | * \param &src matrix to compare to
|
---|
139 | * \return reference to this
|
---|
140 | */
|
---|
141 | MatrixContent &MatrixContent::operator=(const MatrixContent &src)
|
---|
142 | {
|
---|
143 | if(&src!=this){
|
---|
144 | gsl_matrix_memcpy(content,src.content);
|
---|
145 | }
|
---|
146 | return *this;
|
---|
147 | }
|
---|
148 |
|
---|
149 | /** Addition operator.
|
---|
150 | * \param &rhs matrix to add
|
---|
151 | * \return reference to this
|
---|
152 | */
|
---|
153 | const MatrixContent &MatrixContent::operator+=(const MatrixContent &rhs)
|
---|
154 | {
|
---|
155 | gsl_matrix_add(content, rhs.content);
|
---|
156 | return *this;
|
---|
157 | }
|
---|
158 |
|
---|
159 | /** Subtraction operator.
|
---|
160 | * \param &rhs matrix to subtract
|
---|
161 | * \return reference to this
|
---|
162 | */
|
---|
163 | const MatrixContent &MatrixContent::operator-=(const MatrixContent &rhs)
|
---|
164 | {
|
---|
165 | gsl_matrix_sub(content, rhs.content);
|
---|
166 | return *this;
|
---|
167 | }
|
---|
168 |
|
---|
169 | /** Multiplication operator.
|
---|
170 | * Note that here matrix have to have same dimensions.
|
---|
171 | * \param &rhs matrix to multiply with
|
---|
172 | * \return reference to this
|
---|
173 | */
|
---|
174 | const MatrixContent &MatrixContent::operator*=(const MatrixContent &rhs)
|
---|
175 | {
|
---|
176 | ASSERT(rows == rhs.rows,
|
---|
177 | "MatrixContent::operator*=() - row dimension differ: "+toString(rows)+" != "+toString(rhs.rows)+".");
|
---|
178 | ASSERT(columns == rhs.columns,
|
---|
179 | "MatrixContent::operator*=() - columns dimension differ: "+toString(columns)+" != "+toString(rhs.columns)+".");
|
---|
180 | (*this) = (*this)*rhs;
|
---|
181 | return *this;
|
---|
182 | }
|
---|
183 |
|
---|
184 | /** Multiplication with copy operator.
|
---|
185 | * \param &rhs matrix to multiply with
|
---|
186 | * \return reference to newly allocated MatrixContent
|
---|
187 | */
|
---|
188 | const MatrixContent MatrixContent::operator*(const MatrixContent &rhs) const
|
---|
189 | {
|
---|
190 | gsl_matrix *res = gsl_matrix_alloc(rows, columns);
|
---|
191 | gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content, rhs.content, 0.0, res);
|
---|
192 | // gsl_matrix is taken over by constructor, hence no free
|
---|
193 | MatrixContent tmp(res);
|
---|
194 | gsl_matrix_free(res);
|
---|
195 | return tmp;
|
---|
196 | }
|
---|
197 |
|
---|
198 | /** Accessor for manipulating component (i,j).
|
---|
199 | * \param i row number
|
---|
200 | * \param j column number
|
---|
201 | * \return reference to component (i,j)
|
---|
202 | */
|
---|
203 | double &MatrixContent::at(size_t i, size_t j)
|
---|
204 | {
|
---|
205 | ASSERT((i>=0) && (i<rows),
|
---|
206 | "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
|
---|
207 | ASSERT((j>=0) && (j<columns),
|
---|
208 | "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
|
---|
209 | return *gsl_matrix_ptr (content, i, j);
|
---|
210 | }
|
---|
211 |
|
---|
212 | /** Constant accessor for (value of) component (i,j).
|
---|
213 | * \param i row number
|
---|
214 | * \param j column number
|
---|
215 | * \return const component (i,j)
|
---|
216 | */
|
---|
217 | const double MatrixContent::at(size_t i, size_t j) const
|
---|
218 | {
|
---|
219 | ASSERT((i>=0) && (i<rows),
|
---|
220 | "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
|
---|
221 | ASSERT((j>=0) && (j<columns),
|
---|
222 | "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
|
---|
223 | return gsl_matrix_get(content, i, j);
|
---|
224 | }
|
---|
225 |
|
---|
226 | /** Setter for component (i,j).
|
---|
227 | * \param i row numbr
|
---|
228 | * \param j column numnber
|
---|
229 | * \param value value to set componnt (i,j) to
|
---|
230 | */
|
---|
231 | void MatrixContent::set(size_t i, size_t j, const double value)
|
---|
232 | {
|
---|
233 | ASSERT((i>=0) && (i<rows),
|
---|
234 | "MatrixContent::set() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
|
---|
235 | ASSERT((j>=0) && (j<columns),
|
---|
236 | "MatrixContent::set() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
|
---|
237 | gsl_matrix_set(content,i,j,value);
|
---|
238 | }
|
---|
239 |
|
---|
240 | /** Return transposed matrix.
|
---|
241 | * \return new matrix that is transposed of this.
|
---|
242 | */
|
---|
243 | MatrixContent MatrixContent::transpose() const
|
---|
244 | {
|
---|
245 | std::cout << "const MatrixContent::transpose()." << std::endl;
|
---|
246 | gsl_matrix *res = gsl_matrix_alloc(rows, columns);
|
---|
247 | gsl_matrix_transpose_memcpy(res, content);
|
---|
248 | MatrixContent newContent(res);
|
---|
249 | gsl_matrix_free(res);
|
---|
250 | return newContent;
|
---|
251 | }
|
---|
252 |
|
---|
253 | /** Turn this matrix into its transposed.
|
---|
254 | * Note that this is only possible if rows == columns.
|
---|
255 | */
|
---|
256 | MatrixContent &MatrixContent::transpose()
|
---|
257 | {
|
---|
258 | std::cout << "MatrixContent::transpose()." << std::endl;
|
---|
259 | ASSERT( rows == columns,
|
---|
260 | "MatrixContent::transpose() - cannot transpose onto itself as matrix not square: "+toString(rows)+"!="+toString(columns)+"!");
|
---|
261 | double tmp;
|
---|
262 | for (size_t i=0;i<rows;i++)
|
---|
263 | for (size_t j=i+1;j<rows;j++) {
|
---|
264 | tmp = at(j,i);
|
---|
265 | at(j,i) = at(i,j);
|
---|
266 | at(i,j) = tmp;
|
---|
267 | }
|
---|
268 | return *this;
|
---|
269 | }
|
---|
270 |
|
---|
271 | /** Transform the matrix to its eigenbasis ans return resulting eigenvalues.
|
---|
272 | * \warn return vector has to be freed'd
|
---|
273 | * \return gsl_vector pointer to vector of eigenvalues
|
---|
274 | */
|
---|
275 | gsl_vector* MatrixContent::transformToEigenbasis()
|
---|
276 | {
|
---|
277 | ASSERT (rows == columns,
|
---|
278 | "MatrixContent::transformToEigenbasis() - only implemented for square matrices: "+toString(rows)+"!="+toString(columns)+"!");
|
---|
279 | gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(rows);
|
---|
280 | gsl_vector *eval = gsl_vector_alloc(rows);
|
---|
281 | gsl_matrix *evec = gsl_matrix_alloc(rows, rows);
|
---|
282 | gsl_eigen_symmv(content, eval, evec, T);
|
---|
283 | gsl_eigen_symmv_free(T);
|
---|
284 | gsl_matrix_memcpy(content, evec);
|
---|
285 | gsl_matrix_free(evec);
|
---|
286 | return eval;
|
---|
287 | }
|
---|
288 |
|
---|
289 | /** Scalar multiplication operator.
|
---|
290 | * \param factor factor to scale with
|
---|
291 | */
|
---|
292 | const MatrixContent &MatrixContent::operator*=(const double factor)
|
---|
293 | {
|
---|
294 | gsl_matrix_scale(content, factor);
|
---|
295 | return *this;
|
---|
296 | }
|
---|
297 |
|
---|
298 | /** Scalar multiplication and copy operator.
|
---|
299 | * \param factor factor to scale with
|
---|
300 | * \param &mat MatrixContent to scale
|
---|
301 | * \return copied and scaled MatrixContent
|
---|
302 | */
|
---|
303 | const MatrixContent operator*(const double factor,const MatrixContent& mat)
|
---|
304 | {
|
---|
305 | MatrixContent tmp = mat;
|
---|
306 | tmp*=factor;
|
---|
307 | return tmp;
|
---|
308 | }
|
---|
309 |
|
---|
310 | /** Scalar multiplication and copy operator (with operands exchanged).
|
---|
311 | * \param &mat MatrixContent to scale
|
---|
312 | * \param factor factor to scale with
|
---|
313 | * \return copied and scaled MatrixContent
|
---|
314 | */
|
---|
315 | const MatrixContent operator*(const MatrixContent &mat,const double factor)
|
---|
316 | {
|
---|
317 | return factor*mat;
|
---|
318 | }
|
---|
319 |
|
---|
320 | /** Equality operator.
|
---|
321 | * Note that we use numerical sensible checking, i.e. with threshold MYEPSILON.
|
---|
322 | * \param &rhs MatrixContent to checks against
|
---|
323 | */
|
---|
324 | bool MatrixContent::operator==(const MatrixContent &rhs) const
|
---|
325 | {
|
---|
326 | if ((rows == rhs.rows) && (columns == rhs.columns)) {
|
---|
327 | for(int i=rows;i--;){
|
---|
328 | for(int j=columns;j--;){
|
---|
329 | if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
|
---|
330 | return false;
|
---|
331 | }
|
---|
332 | }
|
---|
333 | }
|
---|
334 | return true;
|
---|
335 | }
|
---|
336 | return false;
|
---|
337 | }
|
---|
338 |
|
---|
339 | Vector operator*(const MatrixContent &mat,const Vector &vec)
|
---|
340 | {
|
---|
341 | Vector res;
|
---|
342 | gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, res.content->content);
|
---|
343 | return res;
|
---|
344 | }
|
---|