Annotation of /branches/lamont/src/lib/common/eigen2x2.c

Parent Directory | Revision Log

---

Revision 2081 - (view) (download) (as text)

1 :	jhr	1640	/*! \file eigen2D.c
2 :			*
3 :			* \author Gordon Kindlmann
4 :			*/
5 :
6 :			/*
7 :			* COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :			* All rights reserved.
9 :			*/
10 :
11 :			/*
12 :			Teem: Tools to process and visualize scientific data and images
13 :			Copyright (C) 2011, 2010, 2009, University of Chicago
14 :			Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
15 :			Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
16 :
17 :			This library is free software; you can redistribute it and/or
18 :			modify it under the terms of the GNU Lesser General Public License
19 :			(LGPL) as published by the Free Software Foundation; either
20 :			version 2.1 of the License, or (at your option) any later version.
21 :			The terms of redistributing and/or modifying this software also
22 :			include exceptions to the LGPL that facilitate static linking.
23 :
24 :			This library is distributed in the hope that it will be useful,
25 :			but WITHOUT ANY WARRANTY; without even the implied warranty of
26 :			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 :			Lesser General Public License for more details.
28 :
29 :			You should have received a copy of the GNU Lesser General Public License
30 :			along with this library; if not, write to Free Software Foundation, Inc.,
31 :			51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
32 :			*/
33 :
34 :			#include <Diderot/diderot.h>
35 :
36 :			#define ROOT_DOUBLE 1
37 :			#define ROOT_TWO 2
38 :
39 :			#define SUB(v,i) (((Diderot_union2_t)(v)).r[i])
40 :
41 :			// OpenCL 1.0 does not specify the C representation of the host types
42 :			// for OpenCL vector types (e.g., cl_float4), so we have to handle this
43 :			// mechanism with a macro.
44 :			#if defined(CL_HOST_VECTORS_ARE_STRUCTS)
45 :			# define VSUB(v,i) (v).s[i]
46 :			#elif defined(CL_HOST_VECTORS_ARE_ARRAYS)
47 :			# define VSUB(v,i) (v)[i]
48 :			#else
49 :			# error unable to access elements of host vectors
50 :			#endif
51 :
52 :			/*
53 :			** Eigensolver for symmetric 2x2 matrix:
54 :			**
55 :			** M00 M01
56 :			** M01 M11
57 :			**
58 :			**
59 :			** Return value indicates something about the eigenvalue solution to
60 :			** the quadratic characteristic equation; see ROOT_ #defines above
61 :			**
62 :			** HEY: the numerical precision issues here merit some more scrutiny.
63 :			*/
64 :			int Diderot_evals2x2 (
65 :			Diderot_real_t eval[2],
66 :			const Diderot_real_t _M00, const Diderot_real_t _M01,
67 :			const Diderot_real_t _M11)
68 :			{
69 :			Diderot_real_t mean, Q, D, M00, M01, M11;
70 :			int roots;
71 :
72 :			/* copy the given matrix elements */
73 :			M00 = _M00;
74 :			M01 = _M01;
75 :			M11 = _M11;
76 :
77 :			/*
78 :			** subtract out the eigenvalue mean (will add back to evals later);
79 :			** helps with numerical stability
80 :			*/
81 :			mean = (M00 + M11)/2.0;
82 :			M00 -= mean;
83 :			M11 -= mean;
84 :
85 :			Q = M00 - M11;
86 :			D = 4.0*M01*M01 + Q*Q;
87 :			if (D > EPSILON) {
88 :			/* two distinct roots */
89 :			Diderot_real_t vv;
90 :			vv = SQRT(D)/2.0;
91 :			eval[0] = vv;
92 :			eval[1] = -vv;
93 :			roots = ROOT_TWO;
94 :			}
95 :			else {
96 :			/* double root */
97 :			eval[0] = eval[1] = 0.0;
98 :			roots = ROOT_DOUBLE;
99 :			}
100 :
101 :			/* add back in the eigenvalue mean */
102 :			eval[0] += mean;
103 :			eval[1] += mean;
104 :
105 :			return roots;
106 :			}
107 :
108 :			int Diderot_evecs2x2 (
109 :			Diderot_real_t eval[2], Diderot_vec2_t evec[2],
110 :			const Diderot_real_t _M00, const Diderot_real_t _M01,
111 :			const Diderot_real_t _M11)
112 :			{
113 :			Diderot_real_t mean, Q, D, M00, M01, M11;
114 :			int roots;
115 :
116 :			/* copy the given matrix elements */
117 :			M00 = _M00;
118 :			M01 = _M01;
119 :			M11 = _M11;
120 :
121 :			/*
122 :			** subtract out the eigenvalue mean (will add back to evals later);
123 :			** helps with numerical stability
124 :			*/
125 :			mean = (M00 + M11)/2.0;
126 :			M00 -= mean;
127 :			M11 -= mean;
128 :
129 :			Q = M00 - M11;
130 :			D = 4.0*M01*M01 + Q*Q;
131 :			if (D > EPSILON) {
132 :			/* two distinct roots */
133 :			Diderot_real_t vv;
134 :			Diderot_vec2_t r1, r2;
135 :			vv = SQRT(D)/2.0;
136 :			eval[0] = vv;
137 :			eval[1] = -vv;
138 :			/* null space of T = M - evec[0]*I ==
139 :			[M00 - vv M01 ]
140 :			[ M01 M11 - vv]
141 :			is evec[0], but we know evec[0] and evec[1] are orthogonal,
142 :			so row span of T is evec[1]
143 :			*/
144 :			r1 = vec2(M00 - vv, M01);
145 :			r2 = vec2(M01, M11 - vv);
146 :			if (dot2(r1,r2) > 0.0) {
147 :			evec[1] = vec2(SUB(r1,0)+SUB(r2,0), SUB(r1,1)+SUB(r2,1));
148 :			}
149 :			else {
150 :			evec[1] = vec2(SUB(r1,0)-SUB(r2,0), SUB(r1,1)-SUB(r2,1));
151 :			}
152 :			evec[1] = normalize2(evec[1]);
153 :			evec[0] = vec2(SUB(evec[1],1), -SUB(evec[1],0));
154 :			evec[0] = normalize2(evec[0]);
155 :			roots = ROOT_TWO;
156 :			}
157 :			else {
158 :			/* double root */
159 :			eval[0] = eval[1] = 0.0;
160 :			/* use any basis for eigenvectors */
161 :			evec[0] = vec2(1.0, 0.0);
162 :			evec[1] = vec2(0.0, 1.0);
163 :			roots = ROOT_DOUBLE;
164 :			}
165 :
166 :			/* add back in the eigenvalue mean */
167 :			eval[0] += mean;
168 :			eval[1] += mean;
169 :
170 :			return roots;
171 :			}

---

root@smlnj-gforge.cs.uchicago.edu	ViewVC Help
Powered by ViewVC 1.0.0