1 /*\r
2 * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;\r
3 * "views" mimic the semantic of C++ references: any operation performed\r
4 * on a "view" is actually performed on the "viewed object"\r
5 *\r
6 * Authors:\r
7 * Marco Cecchetti <mrcekets at gmail.com>\r
8 *\r
9 * Copyright 2008 authors\r
10 *\r
11 * This library is free software; you can redistribute it and/or\r
12 * modify it either under the terms of the GNU Lesser General Public\r
13 * License version 2.1 as published by the Free Software Foundation\r
14 * (the "LGPL") or, at your option, under the terms of the Mozilla\r
15 * Public License Version 1.1 (the "MPL"). If you do not alter this\r
16 * notice, a recipient may use your version of this file under either\r
17 * the MPL or the LGPL.\r
18 *\r
19 * You should have received a copy of the LGPL along with this library\r
20 * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
22 * You should have received a copy of the MPL along with this library\r
23 * in the file COPYING-MPL-1.1\r
24 *\r
25 * The contents of this file are subject to the Mozilla Public License\r
26 * Version 1.1 (the "License"); you may not use this file except in\r
27 * compliance with the License. You may obtain a copy of the License at\r
28 * http://www.mozilla.org/MPL/\r
29 *\r
30 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
31 * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
32 * the specific language governing rights and limitations.\r
33 */\r
34 \r
35 \r
36 \r
37 \r
38 #ifndef _NL_MATRIX_H_\r
39 #define _NL_MATRIX_H_\r
40 \r
41 #include <2geom/numeric/vector.h>\r
42 \r
43 #include <cassert>\r
44 #include <utility> // for std::pair\r
45 #include <algorithm> // for std::swap\r
46 #include <sstream>\r
47 #include <string>\r
48 \r
49 #include <gsl/gsl_matrix.h>\r
50 #include <gsl/gsl_linalg.h>\r
51 \r
52 \r
53 namespace Geom { namespace NL {\r
54 \r
55 namespace detail\r
56 {\r
57 \r
58 class BaseMatrixImpl\r
59 {\r
60 public:\r
61 virtual ~BaseMatrixImpl()\r
62 {\r
63 }\r
64 \r
65 ConstVectorView row_const_view(size_t i) const\r
66 {\r
67 return ConstVectorView(gsl_matrix_const_row(m_matrix, i));\r
68 }\r
69 \r
70 ConstVectorView column_const_view(size_t i) const\r
71 {\r
72 return ConstVectorView(gsl_matrix_const_column(m_matrix, i));\r
73 }\r
74 \r
75 const double & operator() (size_t i, size_t j) const\r
76 {\r
77 return *gsl_matrix_const_ptr(m_matrix, i, j);\r
78 }\r
79 \r
80 const gsl_matrix* get_gsl_matrix() const\r
81 {\r
82 return m_matrix;\r
83 }\r
84 \r
85 bool is_zero() const\r
86 {\r
87 return gsl_matrix_isnull(m_matrix);\r
88 }\r
89 \r
90 bool is_positive() const\r
91 {\r
92 return gsl_matrix_ispos(m_matrix);\r
93 }\r
94 \r
95 bool is_negative() const\r
96 {\r
97 return gsl_matrix_isneg(m_matrix);\r
98 }\r
99 \r
100 bool is_non_negative() const\r
101 {\r
102 for ( unsigned int i = 0; i < rows(); ++i )\r
103 {\r
104 for ( unsigned int j = 0; j < columns(); ++j )\r
105 {\r
106 if ( (*this)(i,j) < 0 ) return false;\r
107 }\r
108 }\r
109 return true;\r
110 }\r
111 \r
112 double max() const\r
113 {\r
114 return gsl_matrix_max(m_matrix);\r
115 }\r
116 \r
117 double min() const\r
118 {\r
119 return gsl_matrix_min(m_matrix);\r
120 }\r
121 \r
122 std::pair<size_t, size_t>\r
123 max_index() const\r
124 {\r
125 std::pair<size_t, size_t> indices;\r
126 gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));\r
127 return indices;\r
128 }\r
129 \r
130 std::pair<size_t, size_t>\r
131 min_index() const\r
132 {\r
133 std::pair<size_t, size_t> indices;\r
134 gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));\r
135 return indices;\r
136 }\r
137 \r
138 size_t rows() const\r
139 {\r
140 return m_rows;\r
141 }\r
142 \r
143 size_t columns() const\r
144 {\r
145 return m_columns;\r
146 }\r
147 \r
148 std::string str() const;\r
149 \r
150 protected:\r
151 size_t m_rows, m_columns;\r
152 gsl_matrix* m_matrix;\r
153 \r
154 }; // end class BaseMatrixImpl\r
155 \r
156 \r
157 inline\r
158 bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)\r
159 {\r
160 if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;\r
161 \r
162 for (size_t i = 0; i < m1.rows(); ++i)\r
163 for (size_t j = 0; j < m1.columns(); ++j)\r
164 if (m1(i,j) != m2(i,j)) return false;\r
165 \r
166 return true;\r
167 }\r
168 \r
169 template< class charT >\r
170 inline\r
171 std::basic_ostream<charT> &\r
172 operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)\r
173 {\r
174 if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;\r
175 \r
176 os << "[[" << _matrix(0,0);\r
177 for (size_t j = 1; j < _matrix.columns(); ++j)\r
178 {\r
179 os << ", " << _matrix(0,j);\r
180 }\r
181 os << "]";\r
182 \r
183 for (size_t i = 1; i < _matrix.rows(); ++i)\r
184 {\r
185 os << ", [" << _matrix(i,0);\r
186 for (size_t j = 1; j < _matrix.columns(); ++j)\r
187 {\r
188 os << ", " << _matrix(i,j);\r
189 }\r
190 os << "]";\r
191 }\r
192 os << "]";\r
193 return os;\r
194 }\r
195 \r
196 inline\r
197 std::string BaseMatrixImpl::str() const\r
198 {\r
199 std::ostringstream oss;\r
200 oss << (*this);\r
201 return oss.str();\r
202 }\r
203 \r
204 \r
205 class MatrixImpl : public BaseMatrixImpl\r
206 {\r
207 public:\r
208 \r
209 typedef BaseMatrixImpl base_type;\r
210 \r
211 void set_all( double x )\r
212 {\r
213 gsl_matrix_set_all(m_matrix, x);\r
214 }\r
215 \r
216 void set_identity()\r
217 {\r
218 gsl_matrix_set_identity(m_matrix);\r
219 }\r
220 \r
221 using base_type::operator();\r
222 \r
223 double & operator() (size_t i, size_t j)\r
224 {\r
225 return *gsl_matrix_ptr(m_matrix, i, j);\r
226 }\r
227 \r
228 using base_type::get_gsl_matrix;\r
229 \r
230 gsl_matrix* get_gsl_matrix()\r
231 {\r
232 return m_matrix;\r
233 }\r
234 \r
235 VectorView row_view(size_t i)\r
236 {\r
237 return VectorView(gsl_matrix_row(m_matrix, i));\r
238 }\r
239 \r
240 VectorView column_view(size_t i)\r
241 {\r
242 return VectorView(gsl_matrix_column(m_matrix, i));\r
243 }\r
244 \r
245 void swap_rows(size_t i, size_t j)\r
246 {\r
247 gsl_matrix_swap_rows(m_matrix, i, j);\r
248 }\r
249 \r
250 void swap_columns(size_t i, size_t j)\r
251 {\r
252 gsl_matrix_swap_columns(m_matrix, i, j);\r
253 }\r
254 \r
255 MatrixImpl & transpose()\r
256 {\r
257 assert(columns() == rows());\r
258 gsl_matrix_transpose(m_matrix);\r
259 return (*this);\r
260 }\r
261 \r
262 MatrixImpl & scale(double x)\r
263 {\r
264 gsl_matrix_scale(m_matrix, x);\r
265 return (*this);\r
266 }\r
267 \r
268 MatrixImpl & translate(double x)\r
269 {\r
270 gsl_matrix_add_constant(m_matrix, x);\r
271 return (*this);\r
272 }\r
273 \r
274 MatrixImpl & operator+=(base_type const& _matrix)\r
275 {\r
276 gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());\r
277 return (*this);\r
278 }\r
279 \r
280 MatrixImpl & operator-=(base_type const& _matrix)\r
281 {\r
282 gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());\r
283 return (*this);\r
284 }\r
285 \r
286 }; // end class MatrixImpl\r
287 \r
288 } // end namespace detail\r
289 \r
290 \r
291 using detail::operator==;\r
292 using detail::operator<<;\r
293 \r
294 \r
295 \r
296 \r
297 class Matrix: public detail::MatrixImpl\r
298 {\r
299 public:\r
300 typedef detail::MatrixImpl base_type;\r
301 \r
302 public:\r
303 // the matrix is not inizialized\r
304 Matrix(size_t n1, size_t n2)\r
305 {\r
306 m_rows = n1;\r
307 m_columns = n2;\r
308 m_matrix = gsl_matrix_alloc(n1, n2);\r
309 }\r
310 \r
311 Matrix(size_t n1, size_t n2, double x)\r
312 {\r
313 m_rows = n1;\r
314 m_columns = n2;\r
315 m_matrix = gsl_matrix_alloc(n1, n2);\r
316 gsl_matrix_set_all(m_matrix, x);\r
317 }\r
318 \r
319 Matrix(Matrix const& _matrix)\r
320 : base_type()\r
321 {\r
322 m_rows = _matrix.rows();\r
323 m_columns = _matrix.columns();\r
324 m_matrix = gsl_matrix_alloc(rows(), columns());\r
325 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
326 }\r
327 \r
328 explicit\r
329 Matrix(base_type::base_type const& _matrix)\r
330 {\r
331 m_rows = _matrix.rows();\r
332 m_columns = _matrix.columns();\r
333 m_matrix = gsl_matrix_alloc(rows(), columns());\r
334 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
335 }\r
336 \r
337 Matrix & operator=(Matrix const& _matrix)\r
338 {\r
339 assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
340 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
341 return *this;\r
342 }\r
343 \r
344 Matrix & operator=(base_type::base_type const& _matrix)\r
345 {\r
346 assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
347 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
348 return *this;\r
349 }\r
350 \r
351 virtual ~Matrix()\r
352 {\r
353 gsl_matrix_free(m_matrix);\r
354 }\r
355 \r
356 Matrix & transpose()\r
357 {\r
358 return static_cast<Matrix &>( base_type::transpose() );\r
359 }\r
360 \r
361 Matrix & scale(double x)\r
362 {\r
363 return static_cast<Matrix &>( base_type::scale(x) );\r
364 }\r
365 \r
366 Matrix & translate(double x)\r
367 {\r
368 return static_cast<Matrix &>( base_type::translate(x) );\r
369 }\r
370 \r
371 Matrix & operator+=(base_type::base_type const& _matrix)\r
372 {\r
373 return static_cast<Matrix &>( base_type::operator+=(_matrix) );\r
374 }\r
375 \r
376 Matrix & operator-=(base_type::base_type const& _matrix)\r
377 {\r
378 return static_cast<Matrix &>( base_type::operator-=(_matrix) );\r
379 }\r
380 \r
381 friend\r
382 void swap(Matrix & m1, Matrix & m2);\r
383 friend\r
384 void swap_any(Matrix & m1, Matrix & m2);\r
385 \r
386 }; // end class Matrix\r
387 \r
388 \r
389 // warning! this operation invalidates any view of the passed matrix objects\r
390 inline\r
391 void swap(Matrix & m1, Matrix & m2)\r
392 {\r
393 assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );\r
394 std::swap(m1.m_matrix, m2.m_matrix);\r
395 }\r
396 \r
397 inline\r
398 void swap_any(Matrix & m1, Matrix & m2)\r
399 {\r
400 std::swap(m1.m_matrix, m2.m_matrix);\r
401 std::swap(m1.m_rows, m2.m_rows);\r
402 std::swap(m1.m_columns, m2.m_columns);\r
403 }\r
404 \r
405 \r
406 \r
407 class ConstMatrixView : public detail::BaseMatrixImpl\r
408 {\r
409 public:\r
410 typedef detail::BaseMatrixImpl base_type;\r
411 \r
412 public:\r
413 ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)\r
414 : m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )\r
415 {\r
416 m_rows = n1;\r
417 m_columns = n2;\r
418 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
419 }\r
420 \r
421 ConstMatrixView(const ConstMatrixView & _matrix)\r
422 : base_type(),\r
423 m_matrix_view(_matrix.m_matrix_view)\r
424 {\r
425 m_rows = _matrix.rows();\r
426 m_columns = _matrix.columns();\r
427 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
428 }\r
429 \r
430 ConstMatrixView(const base_type & _matrix)\r
431 : m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))\r
432 {\r
433 m_rows = _matrix.rows();\r
434 m_columns = _matrix.columns();\r
435 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
436 }\r
437 \r
438 private:\r
439 gsl_matrix_const_view m_matrix_view;\r
440 \r
441 }; // end class ConstMatrixView\r
442 \r
443 \r
444 \r
445 \r
446 class MatrixView : public detail::MatrixImpl\r
447 {\r
448 public:\r
449 typedef detail::MatrixImpl base_type;\r
450 \r
451 public:\r
452 MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)\r
453 {\r
454 m_rows = n1;\r
455 m_columns = n2;\r
456 m_matrix_view\r
457 = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);\r
458 m_matrix = &(m_matrix_view.matrix);\r
459 }\r
460 \r
461 MatrixView(const MatrixView & _matrix)\r
462 : base_type()\r
463 {\r
464 m_rows = _matrix.rows();\r
465 m_columns = _matrix.columns();\r
466 m_matrix_view = _matrix.m_matrix_view;\r
467 m_matrix = &(m_matrix_view.matrix);\r
468 }\r
469 \r
470 MatrixView(Matrix & _matrix)\r
471 {\r
472 m_rows = _matrix.rows();\r
473 m_columns = _matrix.columns();\r
474 m_matrix_view\r
475 = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());\r
476 m_matrix = &(m_matrix_view.matrix);\r
477 }\r
478 \r
479 MatrixView & operator=(MatrixView const& _matrix)\r
480 {\r
481 assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
482 gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);\r
483 return *this;\r
484 }\r
485 \r
486 MatrixView & operator=(base_type::base_type const& _matrix)\r
487 {\r
488 assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
489 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
490 return *this;\r
491 }\r
492 \r
493 MatrixView & transpose()\r
494 {\r
495 return static_cast<MatrixView &>( base_type::transpose() );\r
496 }\r
497 \r
498 MatrixView & scale(double x)\r
499 {\r
500 return static_cast<MatrixView &>( base_type::scale(x) );\r
501 }\r
502 \r
503 MatrixView & translate(double x)\r
504 {\r
505 return static_cast<MatrixView &>( base_type::translate(x) );\r
506 }\r
507 \r
508 MatrixView & operator+=(base_type::base_type const& _matrix)\r
509 {\r
510 return static_cast<MatrixView &>( base_type::operator+=(_matrix) );\r
511 }\r
512 \r
513 MatrixView & operator-=(base_type::base_type const& _matrix)\r
514 {\r
515 return static_cast<MatrixView &>( base_type::operator-=(_matrix) );\r
516 }\r
517 \r
518 friend\r
519 void swap_view(MatrixView & m1, MatrixView & m2);\r
520 \r
521 private:\r
522 gsl_matrix_view m_matrix_view;\r
523 \r
524 }; // end class MatrixView\r
525 \r
526 \r
527 inline\r
528 void swap_view(MatrixView & m1, MatrixView & m2)\r
529 {\r
530 assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );\r
531 std::swap(m1.m_matrix_view, m2.m_matrix_view);\r
532 }\r
533 \r
534 Vector operator*( detail::BaseMatrixImpl const& A,\r
535 detail::BaseVectorImpl const& v );\r
536 \r
537 Matrix operator*( detail::BaseMatrixImpl const& A,\r
538 detail::BaseMatrixImpl const& B );\r
539 \r
540 Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);\r
541 \r
542 } } // end namespaces\r
543 \r
544 #endif /*_NL_MATRIX_H_*/\r
545 \r
546 /*\r
547 Local Variables:\r
548 mode:c++\r
549 c-file-style:"stroustrup"\r
550 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
551 indent-tabs-mode:nil\r
552 fill-column:99\r
553 End:\r
554 */\r
555 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r