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

SCM Repository

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

Annotation of /trunk/src/lib/common/eigen2x2.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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

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