Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /branches/cuda/src/compiler/cuda-target/fragments/eigen2x2.in
 [diderot] / branches / cuda / src / compiler / cuda-target / fragments / eigen2x2.in

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

Fri Dec 16 22:17:33 2011 UTC (8 years, 10 months ago) by jhr
File size: 2896 byte(s)
```  Working on CUDA target
```
```#define ROOT_DOUBLE 1
#define ROOT_TWO    2

/*
** Eigensolver for symmetric 2x2 matrix:
**
**  M00  M01
**  M01  M11
**
**
** Return value indicates something about the eigenvalue solution to
** the quadratic characteristic equation; see ROOT_ #defines above
**
** HEY: the numerical precision issues here merit some more scrutiny.
*/
int Diderot_evals2x2 (
float2 * eval,
const Diderot_real_t _M00, const Diderot_real_t _M01,
const Diderot_real_t _M11)
{
Diderot_real_t mean, Q, D, M00, M01, M11;
int roots;

/* copy the given matrix elements */
M00 = _M00;
M01 = _M01;
M11 = _M11;

/*
** subtract out the eigenvalue mean (will add back to evals later);
** helps with numerical stability
*/
mean = (M00 + M11)/2.0;
M00 -= mean;
M11 -= mean;

Q = M00 - M11;
D = 4.0*M01*M01 + Q*Q;
if (D > EPSILON) {
/* two distinct roots */
Diderot_real_t vv;
vv = sqrt(D)/2.0;
eval->s0 = vv;
eval->s1= -vv;
roots = ROOT_TWO;
}
else {
/* double root */
eval->s0 = eval->s1 = 0.0;
roots = ROOT_DOUBLE;
}

/* add back in the eigenvalue mean */
eval->s0 += mean;
eval->s1 += mean;

return roots;
}

int Diderot_evecs2x2 (
float2 * eval, float2 * evec,
const Diderot_real_t _M00, const Diderot_real_t _M01,
const Diderot_real_t _M11)
{
Diderot_real_t mean, Q, D, M00, M01, M11;
int roots;

/* copy the given matrix elements */
M00 = _M00;
M01 = _M01;
M11 = _M11;

/*
** subtract out the eigenvalue mean (will add back to evals later);
** helps with numerical stability
*/
mean = (M00 + M11)/2.0;
M00 -= mean;
M11 -= mean;

Q = M00 - M11;
D = 4.0*M01*M01 + Q*Q;
if (D > EPSILON) {
/* two distinct roots */
Diderot_real_t vv;
Diderot_vec2_t r1, r2;
vv = sqrt(D)/2.0;
eval->s0 = vv;
eval->s1 = -vv;
/* null space of T = M - evec[0]*I ==
[M00 - vv      M01  ]
[  M01      M11 - vv]
is evec[0], but we know evec[0] and evec[1] are orthogonal,
so row span of T is evec[1]
*/
r1 = VEC2(M00 - vv, M01);
r2 = VEC2(M01, M11 - vv);
if (dot(r1,r2) > 0.0) {
evec[1] = VEC2(r1.s0+r2.s0, r1.s1+r2.s1);
}
else {
evec[1] = VEC2(r1.s0-r2.s0, r1.s1-r2.s1);
}
evec[1] = normalize(evec[1]);
evec[0] = VEC2(evec[1].s1, -evec[1].s0);
evec[0] = normalize(evec[0]);
roots = ROOT_TWO;
}
else {
/* double root */
eval->s0 = eval->s1 = 0.0;
/* use any basis for eigenvectors */
evec[0] = VEC2(1.0, 0.0);
evec[1] = VEC2(0.0, 1.0);
roots = ROOT_DOUBLE;
}

/* add back in the eigenvalue mean */
eval->s0 += mean;
eval->s1 += mean;

return roots;
}
```