Bounding Volumes - Real-Time Collision Detection - page 95

done, all rotations are also concatenated into another matrix. Upon exit, this matrix
will contain the eigenvectors. Ideally, this decomposition should be performed in
double-precision arithmetic to minimize numerical errors. The following code for the
Jacobi method is based on the presentation in [Golub96]. First is a subroutine for
assisting in computing the rotation matrix.
// 2-by-2 Symmetric Schur decomposition. Given an n-by-n symmetric matrix
// and indices p, q such that 1 <=p<q<=n,computes a sine-cosine pair
// (s, c) that will serve to form a Jacobi rotation matrix.
//
// See Golub, Van Loan, Matrix Computations, 3rd ed, p428
void SymSchur2(Matrix33 &a, int p, int q, float &c, float &s)
{
if (Abs(a[p][q]) > 0.0001f) {
float r = (a[q][q] - a[p][p]) / (2.0f * a[p][q]);
float t;
if (r >= 0.0f)
t = 1.0f / (r + Sqrt(1.0f + r*r));
else
t = -1.0f / (-r + Sqrt(1.0f + r*r));
c = 1.0f / Sqrt(1.0f + t*t);
s=t*c;
} else {
c = 1.0f;
s = 0.0f;
}
}
Given this support function, the full Jacobi method is now implemented as:
// Computes the eigenvectors and eigenvalues of the symmetric matrix A using
// the classic Jacobi method of iteratively updating A asA=J ∧ T*A*J,
// where J = J(p, q, theta) is the Jacobi rotation matrix.
//
// On exit, v will contain the eigenvectors, and the diagonal elements
// of a are the corresponding eigenvalues.
//
// See Golub, Van Loan, Matrix Computations, 3rd ed, p428
void Jacobi(Matrix33 &a, Matrix33 &v)
{
int i, j, n, p, q;