GCC Code Coverage Report

Directory: ./
File: include/2geom/numeric/symmetric-matrix-fs.h
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 0 52 0.0%
Functions: 0 36 0.0%
Branches: 0 38 0.0%
Line Branch Exec Source
1 /*
2 * SymmetricMatrix, ConstSymmetricMatrixView, SymmetricMatrixView template
3 * classes implement fixed size symmetric matrix; "views" mimic the semantic
4 * of C++ references: any operation performed on a "view" is actually performed
5 * on the "viewed object"
6 *
7 * Authors:
8 * Marco Cecchetti <mrcekets at gmail.com>
9 *
10 * Copyright 2009 authors
11 *
12 * This library is free software; you can redistribute it and/or
13 * modify it either under the terms of the GNU Lesser General Public
14 * License version 2.1 as published by the Free Software Foundation
15 * (the "LGPL") or, at your option, under the terms of the Mozilla
16 * Public License Version 1.1 (the "MPL"). If you do not alter this
17 * notice, a recipient may use your version of this file under either
18 * the MPL or the LGPL.
19 *
20 * You should have received a copy of the LGPL along with this library
21 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 * You should have received a copy of the MPL along with this library
24 * in the file COPYING-MPL-1.1
25 *
26 * The contents of this file are subject to the Mozilla Public License
27 * Version 1.1 (the "License"); you may not use this file except in
28 * compliance with the License. You may obtain a copy of the License at
29 * http://www.mozilla.org/MPL/
30 *
31 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
32 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
33 * the specific language governing rights and limitations.
34 */
35
36
37 #ifndef _NL_SYMMETRIC_MATRIX_FS_H_
38 #define _NL_SYMMETRIC_MATRIX_FS_H_
39
40
41 #include <2geom/numeric/vector.h>
42 #include <2geom/numeric/matrix.h>
43 #include <2geom/utils.h>
44 #include <2geom/exception.h>
45
46 #include <boost/static_assert.hpp>
47
48 #include <cassert>
49 #include <utility> // for std::pair
50 #include <algorithm> // for std::swap, std::copy
51 #include <sstream>
52 #include <string>
53
54
55
56 namespace Geom { namespace NL {
57
58
59 namespace detail
60 {
61
62 template <size_t I, size_t J, bool B = (I < J)>
63 struct index
64 {
65 static const size_t K = index<J, I>::K;
66 };
67
68 template <size_t I, size_t J>
69 struct index<I, J, false>
70 {
71 static const size_t K = (((I+1) * I) >> 1) + J;
72 };
73
74 } // end namespace detail
75
76
77
78
79 template <size_t N>
80 class ConstBaseSymmetricMatrix;
81
82 template <size_t N>
83 class BaseSymmetricMatrix;
84
85 template <size_t N>
86 class SymmetricMatrix;
87
88 template <size_t N>
89 class ConstSymmetricMatrixView;
90
91 template <size_t N>
92 class SymmetricMatrixView;
93
94
95
96 // declaration needed for friend clause
97 template <size_t N>
98 bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1,
99 ConstBaseSymmetricMatrix<N> const& _smatrix2);
100
101
102
103
104 template <size_t N>
105 class ConstBaseSymmetricMatrix
106 {
107 public:
108 const static size_t DIM = N;
109 const static size_t DATA_SIZE = ((DIM+1) * DIM) / 2;
110
111 public:
112
113 ✗ ConstBaseSymmetricMatrix (VectorView const& _data)
114 ✗ : m_data(_data)
115 {
116 ✗ }
117
118 ✗ double operator() (size_t i, size_t j) const
119 {
120 ✗ return m_data[get_index(i,j)];
121 }
122
123 template <size_t I, size_t J>
124 ✗ double get() const
125 {
126 BOOST_STATIC_ASSERT ((I < N && J < N));
127 ✗ return m_data[detail::index<I, J>::K];
128 }
129
130
131 size_t rows() const
132 {
133 return DIM;
134 }
135
136 size_t columns() const
137 {
138 return DIM;
139 }
140
141 ✗ bool is_zero() const
142 {
143 ✗ return m_data.is_zero();
144 }
145
146 bool is_positive() const
147 {
148 return m_data.is_positive();
149 }
150
151 bool is_negative() const
152 {
153 return m_data.is_negative();
154 }
155
156 bool is_non_negative() const
157 {
158 return m_data.is_non_negative();
159 }
160
161 double min() const
162 {
163 return m_data.min();
164 }
165
166 double max() const
167 {
168 return m_data.max();
169 }
170
171 std::pair<size_t, size_t>
172 min_index() const
173 {
174 std::pair<size_t, size_t> indices(0,0);
175 double min_value = m_data[0];
176 for (size_t i = 1; i < DIM; ++i)
177 {
178 for (size_t j = 0; j <= i; ++j)
179 {
180 if (min_value > (*this)(i,j))
181 {
182 min_value = (*this)(i,j);
183 indices.first = i;
184 indices.second = j;
185 }
186 }
187 }
188 return indices;
189 }
190
191 std::pair<size_t, size_t>
192 max_index() const
193 {
194 std::pair<size_t, size_t> indices(0,0);
195 double max_value = m_data[0];
196 for (size_t i = 1; i < DIM; ++i)
197 {
198 for (size_t j = 0; j <= i; ++j)
199 {
200 if (max_value < (*this)(i,j))
201 {
202 max_value = (*this)(i,j);
203 indices.first = i;
204 indices.second = j;
205 }
206 }
207 }
208 return indices;
209 }
210
211 size_t min_on_row_index (size_t i) const
212 {
213 size_t idx = 0;
214 double min_value = (*this)(i,0);
215 for (size_t j = 1; j < DIM; ++j)
216 {
217 if (min_value > (*this)(i,j))
218 {
219 min_value = (*this)(i,j);
220 idx = j;
221 }
222 }
223 return idx;
224 }
225
226 size_t max_on_row_index (size_t i) const
227 {
228 size_t idx = 0;
229 double max_value = (*this)(i,0);
230 for (size_t j = 1; j < DIM; ++j)
231 {
232 if (max_value < (*this)(i,j))
233 {
234 max_value = (*this)(i,j);
235 idx = j;
236 }
237 }
238 return idx;
239 }
240
241 size_t min_on_column_index (size_t j) const
242 {
243 return min_on_row_index(j);
244 }
245
246 size_t max_on_column_index (size_t j) const
247 {
248 return max_on_row_index(j);
249 }
250
251 size_t min_on_diag_index () const
252 {
253 size_t idx = 0;
254 double min_value = (*this)(0,0);
255 for (size_t i = 1; i < DIM; ++i)
256 {
257 if (min_value > (*this)(i,i))
258 {
259 min_value = (*this)(i,i);
260 idx = i;
261 }
262 }
263 return idx;
264 }
265
266 size_t max_on_diag_index () const
267 {
268 size_t idx = 0;
269 double max_value = (*this)(0,0);
270 for (size_t i = 1; i < DIM; ++i)
271 {
272 if (max_value < (*this)(i,i))
273 {
274 max_value = (*this)(i,i);
275 idx = i;
276 }
277 }
278 return idx;
279 }
280
281 std::string str() const;
282
283 ConstSymmetricMatrixView<N-1> main_minor_const_view() const;
284
285 SymmetricMatrix<N> operator- () const;
286
287 ✗ Vector operator* (ConstVectorView _vector) const
288 {
289 ✗ assert (_vector.size() == DIM);
290 ✗ Vector result(DIM, 0.0);
291
292 ✗ for (size_t i = 0; i < DIM; ++i)
293 {
294 ✗ for (size_t j = 0; j < DIM; ++j)
295 {
296 ✗ result[i] += (*this)(i,j) * _vector[j];
297 }
298 }
299 ✗ return result;
300 ✗ }
301
302 protected:
303 static size_t get_index (size_t i, size_t j)
304 {
305 if (i < j) return get_index (j, i);
306 size_t k = (i+1) * i;
307 k >>= 1;
308 k += j;
309 return k;
310 }
311
312 protected:
313 ConstVectorView get_data() const
314 {
315 return m_data;
316 }
317
318 friend
319 bool operator==<N> (ConstBaseSymmetricMatrix const& _smatrix1,
320 ConstBaseSymmetricMatrix const& _smatrix2);
321
322 protected:
323 VectorView m_data;
324
325 }; //end ConstBaseSymmetricMatrix
326
327
328 template <size_t N>
329 class BaseSymmetricMatrix : public ConstBaseSymmetricMatrix<N>
330 {
331 public:
332 typedef ConstBaseSymmetricMatrix<N> base_type;
333
334
335 public:
336
337 ✗ BaseSymmetricMatrix (VectorView const& _data)
338 ✗ : base_type(_data)
339 {
340 ✗ }
341
342 double operator() (size_t i, size_t j) const
343 {
344 return m_data[base_type::get_index(i,j)];
345 }
346
347 ✗ double& operator() (size_t i, size_t j)
348 {
349 ✗ return m_data[base_type::get_index(i,j)];
350 }
351
352 template <size_t I, size_t J>
353 ✗ double& get()
354 {
355 BOOST_STATIC_ASSERT ((I < N && J < N));
356 ✗ return m_data[detail::index<I, J>::K];
357 }
358
359 void set_all (double x)
360 {
361 m_data.set_all(x);
362 }
363
364 SymmetricMatrixView<N-1> main_minor_view();
365
366 BaseSymmetricMatrix& transpose() const
367 {
368 return (*this);
369 }
370
371 BaseSymmetricMatrix& translate (double c)
372 {
373 m_data.translate(c);
374 return (*this);
375 }
376
377 BaseSymmetricMatrix& scale (double c)
378 {
379 m_data.scale(c);
380 return (*this);
381 }
382
383 BaseSymmetricMatrix& operator+= (base_type const& _smatrix)
384 {
385 m_data += (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data);
386 return (*this);
387 }
388
389 BaseSymmetricMatrix& operator-= (base_type const& _smatrix)
390 {
391 m_data -= (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data);
392 return (*this);
393 }
394
395 using base_type::DIM;
396 using base_type::DATA_SIZE;
397 using base_type::m_data;
398 using base_type::operator-;
399 using base_type::operator*;
400
401 }; //end BaseSymmetricMatrix
402
403
404 template <size_t N>
405 class SymmetricMatrix : public BaseSymmetricMatrix<N>
406 {
407 public:
408 typedef BaseSymmetricMatrix<N> base_type;
409 typedef typename base_type::base_type base_base_type;
410
411 using base_type::DIM;
412 using base_type::DATA_SIZE;
413 using base_type::m_data;
414
415 public:
416 ✗ SymmetricMatrix ()
417 ✗ : base_type (VectorView(m_adata, DATA_SIZE))
418 {
419 ✗ }
420
421 explicit
422 SymmetricMatrix (ConstVectorView _data)
423 : base_type (VectorView(m_adata, DATA_SIZE))
424 {
425 assert (_data.size() == DATA_SIZE);
426 m_data = _data;
427 }
428
429 explicit
430 ✗ SymmetricMatrix (const double _data[DATA_SIZE])
431 ✗ : base_type (VectorView(m_adata, DATA_SIZE))
432 {
433 ✗ std::copy (_data, _data + DATA_SIZE, m_adata);
434 ✗ }
435
436 SymmetricMatrix (SymmetricMatrix const& _smatrix)
437 : base_type (VectorView(m_adata, DATA_SIZE))
438 {
439 m_data = _smatrix.m_data;
440 }
441
442 explicit
443 SymmetricMatrix (base_base_type const& _smatrix)
444 : base_type (VectorView(m_adata, DATA_SIZE))
445 {
446 m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
447 }
448
449 explicit
450 SymmetricMatrix (ConstMatrixView const& _matrix)
451 : base_type (VectorView(m_adata, DATA_SIZE))
452 {
453 assert (_matrix.rows() == N && _matrix.columns() == N);
454 for (size_t i = 0; i < N; ++i)
455 for (size_t j = 0; j <= i ; ++j)
456 (*this)(i,j) = _matrix(i,j);
457 }
458
459 SymmetricMatrix& operator= (SymmetricMatrix const& _smatrix)
460 {
461 m_data = _smatrix.m_data;
462 return (*this);
463 }
464
465 SymmetricMatrix& operator= (base_base_type const& _smatrix)
466 {
467
468 m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
469 return (*this);
470 }
471
472 SymmetricMatrix& operator= (ConstMatrixView const& _matrix)
473 {
474 assert (_matrix.rows() == N && _matrix.columns() == N);
475 for (size_t i = 0; i < N; ++i)
476 for (size_t j = 0; j <= i ; ++j)
477 (*this)(i,j) = _matrix(i,j);
478
479 return (*this);
480 }
481
482 // needed for accessing m_adata
483 friend class ConstSymmetricMatrixView<DIM>;
484 friend class SymmetricMatrixView<DIM>;
485 private:
486 double m_adata[DATA_SIZE];
487 }; //end SymmetricMatrix
488
489
490 template <size_t N>
491 class ConstSymmetricMatrixView : public ConstBaseSymmetricMatrix<N>
492 {
493 public:
494 typedef ConstBaseSymmetricMatrix<N> base_type;
495
496 using base_type::DIM;
497 using base_type::DATA_SIZE;
498 using base_type::m_data;
499
500
501 public:
502
503 explicit
504 ✗ ConstSymmetricMatrixView (ConstVectorView _data)
505 ✗ : base_type (const_vector_view_cast(_data))
506 {
507 ✗ assert (_data.size() == DATA_SIZE);
508 ✗ }
509
510 explicit
511 ConstSymmetricMatrixView (const double _data[DATA_SIZE])
512 : base_type (const_vector_view_cast (ConstVectorView (_data, DATA_SIZE)))
513 {
514 }
515
516 ConstSymmetricMatrixView (const ConstSymmetricMatrixView & _smatrix)
517 : base_type (_smatrix.m_data)
518 {
519 }
520
521 ConstSymmetricMatrixView (const base_type & _smatrix)
522 : base_type (static_cast<const ConstSymmetricMatrixView &>(_smatrix).m_data)
523 {
524 }
525
526 }; //end SymmetricMatrix
527
528
529 // declaration needed for friend clause
530 template <size_t N>
531 void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2);
532
533
534 template <size_t N>
535 class SymmetricMatrixView : public BaseSymmetricMatrix<N>
536 {
537 public:
538 typedef BaseSymmetricMatrix<N> base_type;
539 typedef typename base_type::base_type base_base_type;
540
541 using base_type::DIM;
542 using base_type::DATA_SIZE;
543 using base_type::m_data;
544
545 public:
546
547 explicit
548 SymmetricMatrixView (VectorView _data)
549 : base_type (_data)
550 {
551 assert (_data.size() == DATA_SIZE);
552 }
553
554 explicit
555 SymmetricMatrixView (double _data[DATA_SIZE])
556 : base_type (VectorView (_data, DATA_SIZE))
557 {
558 }
559
560 SymmetricMatrixView (const SymmetricMatrixView & _smatrix)
561 : base_type (_smatrix.m_data)
562 {
563 }
564
565 SymmetricMatrixView (SymmetricMatrix<DIM> & _smatrix)
566 : base_type (VectorView (_smatrix.m_adata, DATA_SIZE))
567 {
568 }
569
570 SymmetricMatrixView& operator= (const SymmetricMatrixView & _smatrix)
571 {
572 m_data = _smatrix.m_data;
573 return (*this);
574 }
575
576 SymmetricMatrixView& operator= (const base_base_type & _smatrix)
577 {
578 m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
579 return (*this);
580 }
581
582 friend
583 void swap_view<N>(SymmetricMatrixView & m1, SymmetricMatrixView & m2);
584
585 }; //end SymmetricMatrix
586
587
588
589
590 /*
591 * class ConstBaseSymmetricMatrix methods
592 */
593
594 template <size_t N>
595 inline
596 std::string ConstBaseSymmetricMatrix<N>::str() const
597 {
598 std::ostringstream oss;
599 oss << (*this);
600 return oss.str();
601 }
602
603 template <size_t N>
604 inline
605 ConstSymmetricMatrixView<N-1>
606 ✗ ConstBaseSymmetricMatrix<N>::main_minor_const_view() const
607 {
608 ✗ ConstVectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM);
609 ✗ ConstSymmetricMatrixView<N-1> mm(data);
610 ✗ return mm;
611 ✗ }
612
613 template <size_t N>
614 inline
615 ✗ SymmetricMatrix<N> ConstBaseSymmetricMatrix<N>::operator- () const
616 {
617 ✗ SymmetricMatrix<N> result;
618 ✗ for (size_t i = 0; i < DATA_SIZE; ++i)
619 {
620 ✗ result.m_data[i] = -m_data[i];
621 }
622 ✗ return result;
623 ✗ }
624
625
626 /*
627 * class ConstBaseSymmetricMatrix friend free functions
628 */
629
630 template <size_t N>
631 inline
632 bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1,
633 ConstBaseSymmetricMatrix<N> const& _smatrix2)
634 {
635 return (_smatrix1.m_data == _smatrix2.m_data);
636 }
637
638 /*
639 * class ConstBaseSymmetricMatrix related free functions
640 */
641
642 template< size_t N, class charT >
643 inline
644 std::basic_ostream<charT> &
645 operator<< (std::basic_ostream<charT> & os,
646 const ConstBaseSymmetricMatrix<N> & _matrix)
647 {
648 os << "[[" << _matrix(0,0);
649 for (size_t j = 1; j < N; ++j)
650 {
651 os << ", " << _matrix(0,j);
652 }
653 os << "]";
654 for (size_t i = 1; i < N; ++i)
655 {
656 os << "\n [" << _matrix(i,0);
657 for (size_t j = 1; j < N; ++j)
658 {
659 os << ", " << _matrix(i,j);
660 }
661 os << "]";
662 }
663 os << "]";
664 return os;
665 }
666
667
668 /*
669 * class ConstBaseSymmetricMatrix specialized methods
670 */
671
672 template<>
673 inline
674 ✗ size_t ConstBaseSymmetricMatrix<2>::get_index (size_t i, size_t j)
675 {
676 ✗ return (i+j);
677 }
678
679 template<>
680 inline
681 ✗ size_t ConstBaseSymmetricMatrix<3>::get_index (size_t i, size_t j)
682 {
683 ✗ size_t k = i + j;
684 ✗ if (i == 2 || j == 2) ++k;
685 ✗ return k;
686 }
687
688
689 /*
690 * class BaseSymmetricMatrix methods
691 */
692
693 template <size_t N>
694 inline
695 SymmetricMatrixView<N-1> BaseSymmetricMatrix<N>::main_minor_view()
696 {
697 VectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM);
698 SymmetricMatrixView<N-1> mm(data);
699 return mm;
700 }
701
702
703 /*
704 * class SymmetricMatrixView friend free functions
705 */
706
707 template <size_t N>
708 inline
709 void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2)
710 {
711 swap_view(m1.m_data, m2.m_data);
712 }
713
714 } /* end namespace NL*/ } /* end namespace Geom*/
715
716
717
718
719 #endif // _NL_SYMMETRIC_MATRIX_FS_H_
720
721
722
723
724 /*
725 Local Variables:
726 mode:c++
727 c-file-style:"stroustrup"
728 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
729 indent-tabs-mode:nil
730 fill-column:99
731 End:
732 */
733 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
734