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