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

SCM Repository

[diderot] Annotation of /branches/lamont/src/lib/common/eigen2x2.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1685 - (view) (download) (as text)
Original Path: branches/vis12/src/lib/common/eigen2x2.c

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