in commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/EigenDecomposition.java [563:712]
private void findEigenVectors(final double[][] householderMatrix) {
final double[][]z = householderMatrix.clone();
final int n = main.length;
realEigenvalues = new double[n];
imagEigenvalues = new double[n];
final double[] e = new double[n];
for (int i = 0; i < n - 1; i++) {
realEigenvalues[i] = main[i];
e[i] = secondary[i];
}
realEigenvalues[n - 1] = main[n - 1];
e[n - 1] = 0;
// Determine the largest main and secondary value in absolute term.
double maxAbsoluteValue = 0;
for (int i = 0; i < n; i++) {
if (JdkMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
maxAbsoluteValue = JdkMath.abs(realEigenvalues[i]);
}
if (JdkMath.abs(e[i]) > maxAbsoluteValue) {
maxAbsoluteValue = JdkMath.abs(e[i]);
}
}
// Make null any main and secondary value too small to be significant
if (maxAbsoluteValue != 0) {
for (int i=0; i < n; i++) {
if (JdkMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
realEigenvalues[i] = 0;
}
if (JdkMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
e[i]=0;
}
}
}
for (int j = 0; j < n; j++) {
int its = 0;
int m;
do {
for (m = j; m < n - 1; m++) {
double delta = JdkMath.abs(realEigenvalues[m]) +
JdkMath.abs(realEigenvalues[m + 1]);
if (JdkMath.abs(e[m]) + delta == delta) {
break;
}
}
if (m != j) {
if (its == MAX_ITER) {
throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
MAX_ITER);
}
its++;
double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
double t = JdkMath.sqrt(1 + q * q);
if (q < 0.0) {
q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
} else {
q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
}
double u = 0.0;
double s = 1.0;
double c = 1.0;
int i;
for (i = m - 1; i >= j; i--) {
double p = s * e[i];
double h = c * e[i];
if (JdkMath.abs(p) >= JdkMath.abs(q)) {
c = q / p;
t = JdkMath.sqrt(c * c + 1.0);
e[i + 1] = p * t;
s = 1.0 / t;
c *= s;
} else {
s = p / q;
t = JdkMath.sqrt(s * s + 1.0);
e[i + 1] = q * t;
c = 1.0 / t;
s *= c;
}
if (e[i + 1] == 0.0) {
realEigenvalues[i + 1] -= u;
e[m] = 0.0;
break;
}
q = realEigenvalues[i + 1] - u;
t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
u = s * t;
realEigenvalues[i + 1] = q + u;
q = c * t - h;
for (int ia = 0; ia < n; ia++) {
p = z[ia][i + 1];
z[ia][i + 1] = s * z[ia][i] + c * p;
z[ia][i] = c * z[ia][i] - s * p;
}
}
if (t == 0.0 && i >= j) {
continue;
}
realEigenvalues[j] -= u;
e[j] = q;
e[m] = 0.0;
}
} while (m != j);
}
//Sort the eigen values (and vectors) in increase order
for (int i = 0; i < n; i++) {
int k = i;
double p = realEigenvalues[i];
for (int j = i + 1; j < n; j++) {
if (realEigenvalues[j] > p) {
k = j;
p = realEigenvalues[j];
}
}
if (k != i) {
realEigenvalues[k] = realEigenvalues[i];
realEigenvalues[i] = p;
for (int j = 0; j < n; j++) {
p = z[j][i];
z[j][i] = z[j][k];
z[j][k] = p;
}
}
}
// Determine the largest eigen value in absolute term.
maxAbsoluteValue = 0;
for (int i = 0; i < n; i++) {
if (JdkMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
maxAbsoluteValue=JdkMath.abs(realEigenvalues[i]);
}
}
// Make null any eigen value too small to be significant
if (maxAbsoluteValue != 0.0) {
for (int i=0; i < n; i++) {
if (JdkMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
realEigenvalues[i] = 0;
}
}
}
eigenvectors = new ArrayRealVector[n];
final double[] tmp = new double[n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
tmp[j] = z[j][i];
}
eigenvectors[i] = new ArrayRealVector(tmp);
}
}