GCC Code Coverage Report

Line	Branch	Exec	Source
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			*/
34
35
36
37
38			#ifndef _NL_MATRIX_H_
39			#define _NL_MATRIX_H_
40
41			#include <2geom/exception.h>
42			#include <2geom/numeric/vector.h>
43
44			#include <cassert>
45			#include <utility> // for std::pair
46			#include <algorithm> // for std::swap
47			#include <sstream>
48			#include <string>
49			#include <gsl/gsl_matrix.h>
50			#include <gsl/gsl_linalg.h>
51
52
53			namespace Geom { namespace NL {
54
55			namespace detail
56			{
57
58			class BaseMatrixImpl
59			{
60			public:
61		✗	virtual ~BaseMatrixImpl()
62		✗	{
63		✗	}
64
65			ConstVectorView row_const_view(size_t i) const
66			{
67			return ConstVectorView(gsl_matrix_const_row(m_matrix, i));
68			}
69
70			ConstVectorView column_const_view(size_t i) const
71			{
72			return ConstVectorView(gsl_matrix_const_column(m_matrix, i));
73			}
74
75		✗	const double & operator() (size_t i, size_t j) const
76			{
77		✗	return *gsl_matrix_const_ptr(m_matrix, i, j);
78			}
79
80		✗	const gsl_matrix* get_gsl_matrix() const
81			{
82		✗	return m_matrix;
83			}
84
85			bool is_zero() const
86			{
87			return gsl_matrix_isnull(m_matrix);
88			}
89
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			}
101
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			}
113
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			}
125
126			double max() const
127			{
128			return gsl_matrix_max(m_matrix);
129			}
130
131			double min() const
132			{
133			return gsl_matrix_min(m_matrix);
134			}
135
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			}
143
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			}
151
152		✗	size_t rows() const
153			{
154		✗	return m_rows;
155			}
156
157		✗	size_t columns() const
158			{
159		✗	return m_columns;
160			}
161
162			std::string str() const;
163
164			protected:
165			size_t m_rows, m_columns;
166			gsl_matrix* m_matrix;
167
168			}; // end class BaseMatrixImpl
169
170
171			inline
172			bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)
173			{
174			if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;
175
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;
179
180			return true;
181			}
182
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;
189
190			os << "[[" << _matrix(0,0);
191			for (size_t j = 1; j < _matrix.columns(); ++j)
192			{
193			os << ", " << _matrix(0,j);
194			}
195			os << "]";
196
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			}
209
210			inline
211			std::string BaseMatrixImpl::str() const
212			{
213			std::ostringstream oss;
214			oss << (*this);
215			return oss.str();
216			}
217
218
219			class MatrixImpl : public BaseMatrixImpl
220			{
221			public:
222
223			typedef BaseMatrixImpl base_type;
224
225			void set_all( double x )
226			{
227			gsl_matrix_set_all(m_matrix, x);
228			}
229
230			void set_identity()
231			{
232			gsl_matrix_set_identity(m_matrix);
233			}
234
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			}
240
241		✗	double & operator() (size_t i, size_t j)
242			{
243		✗	return *gsl_matrix_ptr(m_matrix, i, j);
244			}
245
246			using base_type::get_gsl_matrix;
247
248		✗	gsl_matrix* get_gsl_matrix()
249			{
250		✗	return m_matrix;
251			}
252
253		✗	VectorView row_view(size_t i)
254			{
255		✗	return VectorView(gsl_matrix_row(m_matrix, i));
256			}
257
258		✗	VectorView column_view(size_t i)
259			{
260		✗	return VectorView(gsl_matrix_column(m_matrix, i));
261			}
262
263			void swap_rows(size_t i, size_t j)
264			{
265			gsl_matrix_swap_rows(m_matrix, i, j);
266			}
267
268			void swap_columns(size_t i, size_t j)
269			{
270			gsl_matrix_swap_columns(m_matrix, i, j);
271			}
272
273			MatrixImpl & transpose()
274			{
275			assert(columns() == rows());
276			gsl_matrix_transpose(m_matrix);
277			return (*this);
278			}
279
280			MatrixImpl & scale(double x)
281			{
282			gsl_matrix_scale(m_matrix, x);
283			return (*this);
284			}
285
286			MatrixImpl & translate(double x)
287			{
288			gsl_matrix_add_constant(m_matrix, x);
289			return (*this);
290			}
291
292			MatrixImpl & operator+=(base_type const& _matrix)
293			{
294			gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());
295			return (*this);
296			}
297
298			MatrixImpl & operator-=(base_type const& _matrix)
299			{
300			gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());
301			return (*this);
302			}
303
304			}; // end class MatrixImpl
305
306			} // end namespace detail
307
308
309			using detail::operator==;
310			using detail::operator<<;
311
312
313			template <size_t N>
314			class ConstBaseSymmetricMatrix;
315
316
317			class Matrix: public detail::MatrixImpl
318			{
319			public:
320			typedef detail::MatrixImpl base_type;
321
322			public:
323			// the matrix is not initialized
324		✗	Matrix(size_t n1, size_t n2)
325		✗	{
326		✗	m_rows = n1;
327		✗	m_columns = n2;
328		✗	m_matrix = gsl_matrix_alloc(n1, n2);
329		✗	}
330
331		✗	Matrix(size_t n1, size_t n2, double x)
332		✗	{
333		✗	m_rows = n1;
334		✗	m_columns = n2;
335		✗	m_matrix = gsl_matrix_alloc(n1, n2);
336		✗	gsl_matrix_set_all(m_matrix, x);
337		✗	}
338
339			Matrix(Matrix const& _matrix)
340			: base_type()
341			{
342			m_rows = _matrix.rows();
343			m_columns = _matrix.columns();
344			m_matrix = gsl_matrix_alloc(rows(), columns());
345			gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
346			}
347
348			explicit
349		✗	Matrix(base_type::base_type const& _matrix)
350		✗	{
351		✗	m_rows = _matrix.rows();
352		✗	m_columns = _matrix.columns();
353		✗	m_matrix = gsl_matrix_alloc(rows(), columns());
354		✗	gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
355		✗	}
356
357			template <size_t N>
358			explicit
359		✗	Matrix(ConstBaseSymmetricMatrix<N> const& _smatrix)
360		✗	{
361		✗	m_rows = N;
362		✗	m_columns = N;
363		✗	m_matrix = gsl_matrix_alloc(N, N);
364		✗	for (size_t i = 0; i < N; ++i)
365		✗	for (size_t j = 0; j < N ; ++j)
366		✗	(*gsl_matrix_ptr(m_matrix, i, j)) = _smatrix(i,j);
367		✗	}
368
369			Matrix & operator=(Matrix const& _matrix)
370			{
371			assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
372			gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
373			return *this;
374			}
375
376			Matrix & operator=(base_type::base_type const& _matrix)
377			{
378			assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
379			gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
380			return *this;
381			}
382
383			template <size_t N>
384			Matrix & operator=(ConstBaseSymmetricMatrix<N> const& _smatrix)
385			{
386			assert (rows() == N && columns() == N);
387			for (size_t i = 0; i < N; ++i)
388			for (size_t j = 0; j < N ; ++j)
389			(*this)(i,j) = _smatrix(i,j);
390			return *this;
391			}
392
393		✗	~Matrix() override
394		✗	{
395		✗	gsl_matrix_free(m_matrix);
396		✗	}
397
398			Matrix & transpose()
399			{
400			return static_cast<Matrix &>( base_type::transpose() );
401			}
402
403			Matrix & scale(double x)
404			{
405			return static_cast<Matrix &>( base_type::scale(x) );
406			}
407
408			Matrix & translate(double x)
409			{
410			return static_cast<Matrix &>( base_type::translate(x) );
411			}
412
413			Matrix & operator+=(base_type::base_type const& _matrix)
414			{
415			return static_cast<Matrix &>( base_type::operator+=(_matrix) );
416			}
417
418			Matrix & operator-=(base_type::base_type const& _matrix)
419			{
420			return static_cast<Matrix &>( base_type::operator-=(_matrix) );
421			}
422
423			friend
424			void swap(Matrix & m1, Matrix & m2);
425			friend
426			void swap_any(Matrix & m1, Matrix & m2);
427
428			}; // end class Matrix
429
430
431			// warning! this operation invalidates any view of the passed matrix objects
432			inline
433			void swap(Matrix & m1, Matrix & m2)
434			{
435			assert(m1.rows() == m2.rows() && m1.columns() == m2.columns());
436			using std::swap;
437			swap(m1.m_matrix, m2.m_matrix);
438			}
439
440			inline void swap_any(Matrix &m1, Matrix &m2)
441			{
442			using std::swap;
443			swap(m1.m_matrix, m2.m_matrix);
444			swap(m1.m_rows, m2.m_rows);
445			swap(m1.m_columns, m2.m_columns);
446			}
447
448
449
450			class ConstMatrixView : public detail::BaseMatrixImpl
451			{
452			public:
453			typedef detail::BaseMatrixImpl base_type;
454
455			public:
456			ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
457			: m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )
458			{
459			m_rows = n1;
460			m_columns = n2;
461			m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
462			}
463
464			ConstMatrixView(const ConstMatrixView & _matrix)
465			: base_type(),
466			m_matrix_view(_matrix.m_matrix_view)
467			{
468			m_rows = _matrix.rows();
469			m_columns = _matrix.columns();
470			m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
471			}
472
473			ConstMatrixView(const base_type & _matrix)
474			: m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))
475			{
476			m_rows = _matrix.rows();
477			m_columns = _matrix.columns();
478			m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
479			}
480
481			private:
482			gsl_matrix_const_view m_matrix_view;
483
484			}; // end class ConstMatrixView
485
486
487
488
489			class MatrixView : public detail::MatrixImpl
490			{
491			public:
492			typedef detail::MatrixImpl base_type;
493
494			public:
495		✗	MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
496		✗	{
497		✗	m_rows = n1;
498		✗	m_columns = n2;
499			m_matrix_view
500		✗	= gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);
501		✗	m_matrix = &(m_matrix_view.matrix);
502		✗	}
503
504			MatrixView(const MatrixView & _matrix)
505			: base_type()
506			{
507			m_rows = _matrix.rows();
508			m_columns = _matrix.columns();
509			m_matrix_view = _matrix.m_matrix_view;
510			m_matrix = &(m_matrix_view.matrix);
511			}
512
513		✗	MatrixView(Matrix & _matrix)
514		✗	{
515		✗	m_rows = _matrix.rows();
516		✗	m_columns = _matrix.columns();
517			m_matrix_view
518		✗	= gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());
519		✗	m_matrix = &(m_matrix_view.matrix);
520		✗	}
521
522			MatrixView & operator=(MatrixView const& _matrix)
523			{
524			assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
525			gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
526			return *this;
527			}
528
529			MatrixView & operator=(base_type::base_type const& _matrix)
530			{
531			assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
532			gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
533			return *this;
534			}
535
536			MatrixView & transpose()
537			{
538			return static_cast<MatrixView &>( base_type::transpose() );
539			}
540
541			MatrixView & scale(double x)
542			{
543			return static_cast<MatrixView &>( base_type::scale(x) );
544			}
545
546			MatrixView & translate(double x)
547			{
548			return static_cast<MatrixView &>( base_type::translate(x) );
549			}
550
551			MatrixView & operator+=(base_type::base_type const& _matrix)
552			{
553			return static_cast<MatrixView &>( base_type::operator+=(_matrix) );
554			}
555
556			MatrixView & operator-=(base_type::base_type const& _matrix)
557			{
558			return static_cast<MatrixView &>( base_type::operator-=(_matrix) );
559			}
560
561			friend
562			void swap_view(MatrixView & m1, MatrixView & m2);
563
564			private:
565			gsl_matrix_view m_matrix_view;
566
567			}; // end class MatrixView
568
569
570			inline
571			void swap_view(MatrixView & m1, MatrixView & m2)
572			{
573			assert(m1.rows() == m2.rows() && m1.columns() == m2.columns());
574			using std::swap;
575			swap(m1.m_matrix_view, m2.m_matrix_view);
576			}
577
578			Vector operator*( detail::BaseMatrixImpl const& A,
579			detail::BaseVectorImpl const& v );
580
581			Matrix operator*( detail::BaseMatrixImpl const& A,
582			detail::BaseMatrixImpl const& B );
583
584			Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);
585
586			double trace (detail::BaseMatrixImpl const& A);
587
588			double det (detail::BaseMatrixImpl const& A);
589
590			} } // end namespaces
591
592			#endif /*_NL_MATRIX_H_*/
593
594			/*
595			Local Variables:
596			mode:c++
597			c-file-style:"stroustrup"
598			c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
599			indent-tabs-mode:nil
600			fill-column:99
601			End:
602			*/
603			// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
604