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

Parent Directory | 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