GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* SymmetricMatrix trace
3			*
4			* Authors:
5			* Marco Cecchetti <mrcekets at gmail.com>
6			*
7			* Copyright 2009 authors
8			*
9			* This library is free software; you can redistribute it and/or
10			* modify it either under the terms of the GNU Lesser General Public
11			* License version 2.1 as published by the Free Software Foundation
12			* (the "LGPL") or, at your option, under the terms of the Mozilla
13			* Public License Version 1.1 (the "MPL"). If you do not alter this
14			* notice, a recipient may use your version of this file under either
15			* the MPL or the LGPL.
16			*
17			* You should have received a copy of the LGPL along with this library
18			* in the file COPYING-LGPL-2.1; if not, write to the Free Software
19			* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20			* You should have received a copy of the MPL along with this library
21			* in the file COPYING-MPL-1.1
22			*
23			* The contents of this file are subject to the Mozilla Public License
24			* Version 1.1 (the "License"); you may not use this file except in
25			* compliance with the License. You may obtain a copy of the License at
26			* http://www.mozilla.org/MPL/
27			*
28			* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
29			* OF ANY KIND, either express or implied. See the LGPL or the MPL for
30			* the specific language governing rights and limitations.
31			*/
32
33
34			#ifndef _NL_TRACE_H_
35			#define _NL_TRACE_H_
36
37
38			#include <2geom/numeric/matrix.h>
39			#include <2geom/numeric/symmetric-matrix-fs.h>
40
41
42
43
44
45			namespace Geom { namespace NL {
46
47
48			namespace detail
49			{
50
51			/*
52			* helper routines
53			*/
54
55			inline
56		✗	int sgn_prod (int x, int y)
57			{
58		✗	if (x == 0 || y == 0) return 0;
59		✗	if (x == y) return 1;
60		✗	return -1;
61			}
62
63			inline
64		✗	bool abs_less (double x, double y)
65			{
66		✗	return (std::fabs(x) < std::fabs(y));
67			}
68
69
70			/*
71			* trace K-th of symmetric matrix S of order N
72			*/
73			template <size_t K, size_t N>
74			struct trace
75			{
76			static double evaluate(const ConstBaseSymmetricMatrix<N> &S);
77			};
78
79			template <size_t N>
80			struct trace<1,N>
81			{
82			static
83			double evaluate (const ConstBaseSymmetricMatrix<N> & S)
84			{
85			double t = 0;
86			for (size_t i = 0; i < N; ++i)
87			{
88			t += S(i,i);
89			}
90			return t;
91			}
92			};
93
94			template <size_t N>
95			struct trace<N,N>
96			{
97			static
98			double evaluate (const ConstBaseSymmetricMatrix<N> & S)
99			{
100			Matrix M(S);
101			return det(M);
102			}
103			};
104
105			/*
106			* trace for symmetric matrix of order 2
107			*/
108			template <>
109			struct trace<1,2>
110			{
111			static
112		✗	double evaluate (const ConstBaseSymmetricMatrix<2> & S)
113			{
114		✗	return (S.get<0,0>() + S.get<1,1>());
115			}
116			};
117
118			template <>
119			struct trace<2,2>
120			{
121			static
122		✗	double evaluate (const ConstBaseSymmetricMatrix<2> & S)
123			{
124		✗	return (S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>());
125			}
126			};
127
128
129			/*
130			* trace for symmetric matrix of order 3
131			*/
132			template <>
133			struct trace<1,3>
134			{
135			static
136			double evaluate (const ConstBaseSymmetricMatrix<3> & S)
137			{
138			return (S.get<0,0>() + S.get<1,1>() + S.get<2,2>());
139			}
140			};
141
142			template <>
143			struct trace<2,3>
144			{
145			static
146			double evaluate (const ConstBaseSymmetricMatrix<3> & S)
147			{
148			double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
149			double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
150			double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
151			return (a00 + a11 + a22);
152			}
153			};
154
155			template <>
156			struct trace<3,3>
157			{
158			static
159			double evaluate (const ConstBaseSymmetricMatrix<3> & S)
160			{
161			double d = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
162			d += (2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>());
163			d -= (S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
164			d -= (S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
165			d -= (S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
166			return d;
167			}
168			};
169
170
171			/*
172			* sign of trace K-th
173			*/
174			template <size_t K, size_t N>
175			struct trace_sgn
176			{
177			static
178		✗	int evaluate (const ConstBaseSymmetricMatrix<N> & S)
179			{
180		✗	double d = trace<K, N>::evaluate(S);
181		✗	return sgn(d);
182			}
183			};
184
185
186			/*
187			* sign of trace for symmetric matrix of order 2
188			*/
189			template <>
190			struct trace_sgn<2,2>
191			{
192			static
193		✗	int evaluate (const ConstBaseSymmetricMatrix<2> & S)
194			{
195		✗	double m00 = S.get<0,0>();
196		✗	double m10 = S.get<1,0>();
197		✗	double m11 = S.get<1,1>();
198
199		✗	int sm00 = sgn (m00);
200		✗	int sm10 = sgn (m10);
201		✗	int sm11 = sgn (m11);
202
203		✗	if (sm10 == 0)
204			{
205		✗	return sgn_prod (sm00, sm11);
206			}
207			else
208			{
209		✗	int sm00m11 = sgn_prod (sm00, sm11);
210		✗	if (sm00m11 == 1)
211			{
212		✗	int e00, e10, e11;
213		✗	double f00 = std::frexp (m00, &e00);
214		✗	double f10 = std::frexp (m10, &e10);
215		✗	double f11 = std::frexp (m11, &e11);
216
217		✗	int e0011 = e00 + e11;
218		✗	int e1010 = e10 << 1;
219		✗	int ed = e0011 - e1010;
220
221		✗	if (ed > 1)
222			{
223		✗	return 1;
224			}
225		✗	else if (ed < -1)
226			{
227		✗	return -1;
228			}
229			else
230			{
231		✗	double d = std::ldexp (f00 * f11, ed) - f10 * f10;
232			//std::cout << "trace_sgn<2,2>: det = " << d << std::endl;
233		✗	double eps = std::ldexp (1, -50);
234		✗	if (std::fabs(d) < eps) return 0;
235		✗	return sgn (d);
236			}
237			}
238		✗	return -1;
239			}
240			}
241			};
242
243
244			/*
245			* sign of trace for symmetric matrix of order 3
246			*/
247			template <>
248			struct trace_sgn<2,3>
249			{
250			static
251		✗	int evaluate (const ConstBaseSymmetricMatrix<3> & S)
252			{
253		✗	double eps = std::ldexp (1, -50);
254		✗	double t[6];
255
256		✗	t[0] = S.get<1,1>() * S.get<2,2>();
257		✗	t[1] = - S.get<1,2>() * S.get<2,1>();
258		✗	t[2] = S.get<0,0>() * S.get<2,2>();
259		✗	t[3] = - S.get<0,2>() * S.get<2,0>();
260		✗	t[4] = S.get<0,0>() * S.get<1,1>();
261		✗	t[5] = - S.get<0,1>() * S.get<1,0>();
262
263
264		✗	double* maxp = std::max_element (t, t+6, abs_less);
265		✗	int em;
266		✗	std::frexp(*maxp, &em);
267		✗	double d = 0;
268		✗	for (double i : t)
269			{
270		✗	d += i;
271			}
272		✗	double r = std::fabs (std::ldexp (d, -em)); // relative error
273			//std::cout << "trace_sgn<2,3>: d = " << d << std::endl;
274			//std::cout << "trace_sgn<2,3>: r = " << r << std::endl;
275		✗	if (r < eps) return 0;
276		✗	if (d > 0) return 1;
277		✗	return -1;
278			}
279			};
280
281			template <>
282			struct trace_sgn<3,3>
283			{
284			static
285		✗	int evaluate (const ConstBaseSymmetricMatrix<3> & S)
286			{
287
288		✗	double eps = std::ldexp (1, -48);
289		✗	double t[5];
290
291		✗	t[0] = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
292		✗	t[1] = 2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>();
293		✗	t[2] = -(S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
294		✗	t[3] = -(S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
295		✗	t[4] = -(S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
296
297		✗	double* maxp = std::max_element (t, t+5, abs_less);
298		✗	int em;
299		✗	std::frexp(*maxp, &em);
300		✗	double d = 0;
301		✗	for (double i : t)
302			{
303		✗	d += i;
304			}
305			//std::cout << "trace_sgn<3,3>: d = " << d << std::endl;
306		✗	double r = std::fabs (std::ldexp (d, -em)); // relative error
307			//std::cout << "trace_sgn<3,3>: r = " << r << std::endl;
308
309		✗	if (r < eps) return 0;
310		✗	if (d > 0) return 1;
311		✗	return -1;
312			}
313			}; // end struct trace_sgn<3,3>
314
315			} // end namespace detail
316
317
318			template <size_t K, size_t N>
319			inline
320			double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
321			{
322			return detail::trace<K, N>::evaluate(_matrix);
323			}
324
325			template <size_t N>
326			inline
327		✗	double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
328			{
329		✗	return detail::trace<1, N>::evaluate(_matrix);
330			}
331
332			template <size_t N>
333			inline
334		✗	double det (const ConstBaseSymmetricMatrix<N> & _matrix)
335			{
336		✗	return detail::trace<N, N>::evaluate(_matrix);
337			}
338
339
340			template <size_t K, size_t N>
341			inline
342		✗	int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
343			{
344		✗	return detail::trace_sgn<K, N>::evaluate(_matrix);
345			}
346
347			template <size_t N>
348			inline
349		✗	int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
350			{
351		✗	return detail::trace_sgn<1, N>::evaluate(_matrix);
352			}
353
354			template <size_t N>
355			inline
356		✗	int det_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
357			{
358		✗	return detail::trace_sgn<N, N>::evaluate(_matrix);
359			}
360
361			/*
362			template <size_t N>
363			inline
364			size_t rank (const ConstBaseSymmetricMatrix<N> & S)
365			{
366			THROW_NOTIMPLEMENTED();
367			return 0;
368			}
369
370			template <>
371			inline
372			size_t rank<2> (const ConstBaseSymmetricMatrix<2> & S)
373			{
374			if (S.is_zero()) return 0;
375			double d = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
376			if (d != 0) return 2;
377			return 1;
378			}
379
380			template <>
381			inline
382			size_t rank<3> (const ConstBaseSymmetricMatrix<3> & S)
383			{
384			if (S.is_zero()) return 0;
385
386			double a20 = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>();
387			double a21 = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>();
388			double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
389			double d = a20 * S.get<2,0>() + a21 * S.get<2,1>() + a22 * S.get<2,2>();
390
391			if (d != 0) return 3;
392
393			if (a20 != 0 || a21 != 0 || a22 != 0) return 2;
394
395			double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
396			if (a00 != 0) return 2;
397
398			double a10 = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>();
399			if (a10 != 0) return 2;
400
401			double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
402			if (a11 != 0) return 2;
403
404			return 1;
405			}
406			*/
407
408			} /* end namespace NL*/ } /* end namespace Geom*/
409
410
411
412
413			#endif // _NL_TRACE_H_
414
415
416
417
418			/*
419			Local Variables:
420			mode:c++
421			c-file-style:"stroustrup"
422			c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
423			indent-tabs-mode:nil
424			fill-column:99
425			End:
426			*/
427			// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
428