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

SCM Repository

[diderot] Annotation of /branches/cuda/src/compiler/cuda-target/fragments/eigen2x2.in
ViewVC logotype

Annotation of /branches/cuda/src/compiler/cuda-target/fragments/eigen2x2.in

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1677 - (view) (download)

1 : jhr 1677 #define ROOT_DOUBLE 1
2 :     #define ROOT_TWO 2
3 :    
4 :     /*
5 :     ** Eigensolver for symmetric 2x2 matrix:
6 :     **
7 :     ** M00 M01
8 :     ** M01 M11
9 :     **
10 :     **
11 :     ** Return value indicates something about the eigenvalue solution to
12 :     ** the quadratic characteristic equation; see ROOT_ #defines above
13 :     **
14 :     ** HEY: the numerical precision issues here merit some more scrutiny.
15 :     */
16 :     int Diderot_evals2x2 (
17 :     float2 * eval,
18 :     const Diderot_real_t _M00, const Diderot_real_t _M01,
19 :     const Diderot_real_t _M11)
20 :     {
21 :     Diderot_real_t mean, Q, D, M00, M01, M11;
22 :     int roots;
23 :    
24 :     /* copy the given matrix elements */
25 :     M00 = _M00;
26 :     M01 = _M01;
27 :     M11 = _M11;
28 :    
29 :     /*
30 :     ** subtract out the eigenvalue mean (will add back to evals later);
31 :     ** helps with numerical stability
32 :     */
33 :     mean = (M00 + M11)/2.0;
34 :     M00 -= mean;
35 :     M11 -= mean;
36 :    
37 :     Q = M00 - M11;
38 :     D = 4.0*M01*M01 + Q*Q;
39 :     if (D > EPSILON) {
40 :     /* two distinct roots */
41 :     Diderot_real_t vv;
42 :     vv = sqrt(D)/2.0;
43 :     eval->s0 = vv;
44 :     eval->s1= -vv;
45 :     roots = ROOT_TWO;
46 :     }
47 :     else {
48 :     /* double root */
49 :     eval->s0 = eval->s1 = 0.0;
50 :     roots = ROOT_DOUBLE;
51 :     }
52 :    
53 :     /* add back in the eigenvalue mean */
54 :     eval->s0 += mean;
55 :     eval->s1 += mean;
56 :    
57 :     return roots;
58 :     }
59 :    
60 :     int Diderot_evecs2x2 (
61 :     float2 * eval, float2 * evec,
62 :     const Diderot_real_t _M00, const Diderot_real_t _M01,
63 :     const Diderot_real_t _M11)
64 :     {
65 :     Diderot_real_t mean, Q, D, M00, M01, M11;
66 :     int roots;
67 :    
68 :     /* copy the given matrix elements */
69 :     M00 = _M00;
70 :     M01 = _M01;
71 :     M11 = _M11;
72 :    
73 :     /*
74 :     ** subtract out the eigenvalue mean (will add back to evals later);
75 :     ** helps with numerical stability
76 :     */
77 :     mean = (M00 + M11)/2.0;
78 :     M00 -= mean;
79 :     M11 -= mean;
80 :    
81 :     Q = M00 - M11;
82 :     D = 4.0*M01*M01 + Q*Q;
83 :     if (D > EPSILON) {
84 :     /* two distinct roots */
85 :     Diderot_real_t vv;
86 :     Diderot_vec2_t r1, r2;
87 :     vv = sqrt(D)/2.0;
88 :     eval->s0 = vv;
89 :     eval->s1 = -vv;
90 :     /* null space of T = M - evec[0]*I ==
91 :     [M00 - vv M01 ]
92 :     [ M01 M11 - vv]
93 :     is evec[0], but we know evec[0] and evec[1] are orthogonal,
94 :     so row span of T is evec[1]
95 :     */
96 :     r1 = VEC2(M00 - vv, M01);
97 :     r2 = VEC2(M01, M11 - vv);
98 :     if (dot(r1,r2) > 0.0) {
99 :     evec[1] = VEC2(r1.s0+r2.s0, r1.s1+r2.s1);
100 :     }
101 :     else {
102 :     evec[1] = VEC2(r1.s0-r2.s0, r1.s1-r2.s1);
103 :     }
104 :     evec[1] = normalize(evec[1]);
105 :     evec[0] = VEC2(evec[1].s1, -evec[1].s0);
106 :     evec[0] = normalize(evec[0]);
107 :     roots = ROOT_TWO;
108 :     }
109 :     else {
110 :     /* double root */
111 :     eval->s0 = eval->s1 = 0.0;
112 :     /* use any basis for eigenvectors */
113 :     evec[0] = VEC2(1.0, 0.0);
114 :     evec[1] = VEC2(0.0, 1.0);
115 :     roots = ROOT_DOUBLE;
116 :     }
117 :    
118 :     /* add back in the eigenvalue mean */
119 :     eval->s0 += mean;
120 :     eval->s1 += mean;
121 :    
122 :     return roots;
123 :     }

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