Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /branches/pure-cfg/src/lib/c-target/eigen2x2.c
ViewVC logotype

Annotation of /branches/pure-cfg/src/lib/c-target/eigen2x2.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : glk 1540 /*! \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 :     #ifdef DIDEROT_DOUBLE_PRECISION
40 : jhr 1544 # define EPSILON 1.0E-12
41 :     # define SQRT(x) sqrt(x)
42 :     # define VEC(a,b) vec2d(a,b)
43 :     # define DOT(u,v) dot2d(u,v)
44 :     # define SCALE(s,v) scale2d(s,v)
45 :     # define NORMALIZE(v) normalize2d(v)
46 : glk 1540 #else
47 : jhr 1544 # define EPSILON 1.0E-12 // FIXME: probably want a different value of epsilon for single-precision
48 :     # define SQRT(x) sqrtf(x)
49 :     # define VEC(a,b) vec2f(a,b)
50 :     # define DOT(u,v) dot2f(u,v)
51 :     # define SCALE(s,v) scale2f(s,v)
52 :     # define NORMALIZE(v) normalize2f(v)
53 : glk 1540 #endif
54 : jhr 1544 #define SUB(v,i) (((Diderot_union2_t)(v)).r[i])
55 : glk 1540
56 :     // OpenCL 1.0 does not specify the C representation of the host types
57 :     // for OpenCL vector types (e.g., cl_float4), so we have to handle this
58 :     // mechanism with a macro.
59 :     #if defined(CL_HOST_VECTORS_ARE_STRUCTS)
60 :     # define VSUB(v,i) (v).s[i]
61 :     #elif defined(CL_HOST_VECTORS_ARE_ARRAYS)
62 :     # define VSUB(v,i) (v)[i]
63 :     #else
64 :     # error unable to access elements of host vectors
65 :     #endif
66 :    
67 :     /*
68 :     ** Eigensolver for symmetric 2x2 matrix:
69 :     **
70 :     ** M00 M01
71 :     ** M01 M11
72 :     **
73 :     **
74 :     ** Return value indicates something about the eigenvalue solution to
75 :     ** the quadratic characteristic equation; see ROOT_ #defines above
76 :     **
77 :     ** HEY: the numerical precision issues here merit some more scrutiny.
78 :     */
79 : jhr 1544 int Diderot_evals2x2 (
80 :     Diderot_real_t eval[2],
81 :     const Diderot_real_t _M00, const Diderot_real_t _M01,
82 :     const Diderot_real_t _M11)
83 : glk 1540 {
84 :     Diderot_real_t mean, Q, D, M00, M01, M11;
85 :     int roots;
86 :    
87 :     /* copy the given matrix elements */
88 :     M00 = _M00;
89 :     M01 = _M01;
90 :     M11 = _M11;
91 :    
92 :     /*
93 :     ** subtract out the eigenvalue mean (will add back to evals later);
94 :     ** helps with numerical stability
95 :     */
96 :     mean = (M00 + M11)/2.0;
97 :     M00 -= mean;
98 :     M11 -= mean;
99 :    
100 :     Q = M00 - M11;
101 :     D = 4.0*M01*M01 + Q*Q;
102 : jhr 1544 if (D > EPSILON) {
103 : glk 1540 /* two distinct roots */
104 : jhr 1544 Diderot_real_t vv;
105 :     vv = SQRT(D)/2.0;
106 :     eval[0] = vv;
107 :     eval[1] = -vv;
108 :     roots = ROOT_TWO;
109 :     }
110 :     else {
111 : glk 1540 /* double root */
112 :     eval[0] = eval[1] = 0.0;
113 :     roots = ROOT_DOUBLE;
114 :     }
115 :    
116 :     /* add back in the eigenvalue mean */
117 :     eval[0] += mean;
118 :     eval[1] += mean;
119 :    
120 :     return roots;
121 :     }
122 :    
123 : jhr 1544 int Diderot_evecs2x2 (
124 :     Diderot_real_t eval[2], Diderot_vec2_t evec[2],
125 :     const Diderot_real_t _M00, const Diderot_real_t _M01,
126 :     const Diderot_real_t _M11)
127 : glk 1540 {
128 :     Diderot_real_t mean, Q, D, M00, M01, M11;
129 :     int roots;
130 :    
131 :     /* copy the given matrix elements */
132 :     M00 = _M00;
133 :     M01 = _M01;
134 :     M11 = _M11;
135 :    
136 :     /*
137 :     ** subtract out the eigenvalue mean (will add back to evals later);
138 :     ** helps with numerical stability
139 :     */
140 :     mean = (M00 + M11)/2.0;
141 :     M00 -= mean;
142 :     M11 -= mean;
143 :    
144 :     Q = M00 - M11;
145 :     D = 4.0*M01*M01 + Q*Q;
146 : jhr 1544 if (D > EPSILON) {
147 : glk 1540 /* two distinct roots */
148 : jhr 1544 Diderot_real_t vv;
149 :     Diderot_vec2_t r1, r2;
150 :     vv = SQRT(D)/2.0;
151 :     eval[0] = vv;
152 :     eval[1] = -vv;
153 : glk 1540 /* null space of T = M - evec[0]*I ==
154 :     [M00 - vv M01 ]
155 :     [ M01 M11 - vv]
156 :     is evec[0], but we know evec[0] and evec[1] are orthogonal,
157 :     so row span of T is evec[1]
158 :     */
159 : jhr 1544 r1 = VEC(M00 - vv, M01);
160 :     r2 = VEC(M01, M11 - vv);
161 :     if (DOT(r1,r2) > 0.0) {
162 :     evec[1] = VEC(SUB(r1,0)+SUB(r2,0), SUB(r1,1)+SUB(r2,1));
163 :     }
164 :     else {
165 :     evec[1] = VEC(SUB(r1,0)-SUB(r2,0), SUB(r1,1)-SUB(r2,1));
166 :     }
167 :     evec[1] = NORMALIZE(evec[1]);
168 :     evec[0] = VEC(SUB(evec[1],1), -SUB(evec[1],0));
169 :     evec[0] = NORMALIZE(evec[0]);
170 :     roots = ROOT_TWO;
171 :     }
172 :     else {
173 : glk 1540 /* double root */
174 : jhr 1544 eval[0] = eval[1] = 0.0;
175 : glk 1540 /* use any basis for eigenvectors */
176 : jhr 1544 evec[0] = VEC(1.0, 0.0);
177 :     evec[1] = VEC(0.0, 1.0);
178 :     roots = ROOT_DOUBLE;
179 : glk 1540 }
180 :    
181 :     /* add back in the eigenvalue mean */
182 :     eval[0] += mean;
183 :     eval[1] += mean;
184 :    
185 :     return roots;
186 :     }

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