Game Development Reference
In-Depth Information
inline __host__ __device__
void solve_eig_psd( float E, float F, float G, float & lambda1,
float & lambda2, float2 &ev)
{
float B=(E+G)/2;
if (B > 0) {
float D=(E-G)/2;
float FF = F*F;
float R= sqrtf (D*D + FF);
lambda1 = B + R;
lambda2 = fmaxf (0, E*G - FF) / lambda1;
if (R>0){
if (D >= 0) {
float nx=D+R;
ev = make_float2 (nx, F) * rsqrtf (nx*nx + FF);
} else {
float ny=-D+R;
ev = make_float2 (F, ny) * rsqrtf (FF + ny*ny);
}
} else {
ev = make_float2 (1, 0);
}
} else {
lambda1 = lambda2 = 0;
ev = make_float2 (1, 0);
}
}
Listing 5.1. Eigenanalysis of a symmetric positive semidefinite 2 × 2matrix.
in loss of accuracy due to cancelation. In case of the major eigenvalue, only
the subtractions under the square root are an issue, because of the assumption
that E + G> 0.
G ) 2 ,since
then catastrophic cancelation is avoided and replaced with benign cancelation
[Goldberg 91]:
The subtractions are best implemented as ( E
E
2 + F 2 .
E + G
2
G
λ 1
=
+ R,
R =
2
For the computation of the minor eigenvalue, we have to take care of the sub-
traction in front of the square root. A common approach is to use the property
that the product of the eigenvalues equals its determinant [Goldberg 91]:
F 2
det( J )
λ 1
EG
λ 2 =
=
.
λ 1
Notice that λ 1
> 0, since we assumed E + G> 0.
Because of rounding er-
F 2 may become negative, and therefore it is advisable to compute
rors, EG
F 2 ).
max(0 ,EG