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 1985 - (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 :     /*
42 :     ** Eigensolver for symmetric 2x2 matrix:
43 :     **
44 :     ** M00 M01
45 :     ** M01 M11
46 :     **
47 :     **
48 :     ** Return value indicates something about the eigenvalue solution to
49 :     ** the quadratic characteristic equation; see ROOT_ #defines above
50 :     **
51 :     ** HEY: the numerical precision issues here merit some more scrutiny.
52 :     */
53 :     int Diderot_evals2x2 (
54 :     Diderot_real_t eval[2],
55 :     const Diderot_real_t _M00, const Diderot_real_t _M01,
56 :     const Diderot_real_t _M11)
57 :     {
58 :     Diderot_real_t mean, Q, D, M00, M01, M11;
59 :     int roots;
60 :    
61 :     /* copy the given matrix elements */
62 :     M00 = _M00;
63 :     M01 = _M01;
64 :     M11 = _M11;
65 :    
66 :     /*
67 :     ** subtract out the eigenvalue mean (will add back to evals later);
68 :     ** helps with numerical stability
69 :     */
70 :     mean = (M00 + M11)/2.0;
71 :     M00 -= mean;
72 :     M11 -= mean;
73 :    
74 :     Q = M00 - M11;
75 :     D = 4.0*M01*M01 + Q*Q;
76 :     if (D > EPSILON) {
77 :     /* two distinct roots */
78 :     Diderot_real_t vv;
79 :     vv = SQRT(D)/2.0;
80 :     eval[0] = vv;
81 :     eval[1] = -vv;
82 :     roots = ROOT_TWO;
83 :     }
84 :     else {
85 :     /* double root */
86 :     eval[0] = eval[1] = 0.0;
87 :     roots = ROOT_DOUBLE;
88 :     }
89 :    
90 :     /* add back in the eigenvalue mean */
91 :     eval[0] += mean;
92 :     eval[1] += mean;
93 :    
94 :     return roots;
95 :     }
96 :    
97 :     int Diderot_evecs2x2 (
98 :     Diderot_real_t eval[2], Diderot_vec2_t evec[2],
99 :     const Diderot_real_t _M00, const Diderot_real_t _M01,
100 :     const Diderot_real_t _M11)
101 :     {
102 :     Diderot_real_t mean, Q, D, M00, M01, M11;
103 :     int roots;
104 :    
105 :     /* copy the given matrix elements */
106 :     M00 = _M00;
107 :     M01 = _M01;
108 :     M11 = _M11;
109 :    
110 :     /*
111 :     ** subtract out the eigenvalue mean (will add back to evals later);
112 :     ** helps with numerical stability
113 :     */
114 :     mean = (M00 + M11)/2.0;
115 :     M00 -= mean;
116 :     M11 -= mean;
117 :    
118 :     Q = M00 - M11;
119 :     D = 4.0*M01*M01 + Q*Q;
120 :     if (D > EPSILON) {
121 :     /* two distinct roots */
122 :     Diderot_real_t vv;
123 :     Diderot_vec2_t r1, r2;
124 :     vv = SQRT(D)/2.0;
125 :     eval[0] = vv;
126 :     eval[1] = -vv;
127 :     /* null space of T = M - evec[0]*I ==
128 :     [M00 - vv M01 ]
129 :     [ M01 M11 - vv]
130 :     is evec[0], but we know evec[0] and evec[1] are orthogonal,
131 :     so row span of T is evec[1]
132 :     */
133 :     r1 = vec2(M00 - vv, M01);
134 :     r2 = vec2(M01, M11 - vv);
135 :     if (dot2(r1,r2) > 0.0) {
136 :     evec[1] = vec2(SUB(r1,0)+SUB(r2,0), SUB(r1,1)+SUB(r2,1));
137 :     }
138 :     else {
139 :     evec[1] = vec2(SUB(r1,0)-SUB(r2,0), SUB(r1,1)-SUB(r2,1));
140 :     }
141 :     evec[1] = normalize2(evec[1]);
142 :     evec[0] = vec2(SUB(evec[1],1), -SUB(evec[1],0));
143 :     evec[0] = normalize2(evec[0]);
144 :     roots = ROOT_TWO;
145 :     }
146 :     else {
147 :     /* double root */
148 :     eval[0] = eval[1] = 0.0;
149 :     /* use any basis for eigenvectors */
150 :     evec[0] = vec2(1.0, 0.0);
151 :     evec[1] = vec2(0.0, 1.0);
152 :     roots = ROOT_DOUBLE;
153 :     }
154 :    
155 :     /* add back in the eigenvalue mean */
156 :     eval[0] += mean;
157 :     eval[1] += mean;
158 :    
159 :     return roots;
160 :     }

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