in commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/EigenDecomposition.java [765:940]
private void findEigenVectorsFromSchur(final SchurTransformer schur)
throws MathArithmeticException {
final double[][] matrixT = schur.getT().getData();
final double[][] matrixP = schur.getP().getData();
final int n = matrixT.length;
// compute matrix norm
double norm = 0.0;
for (int i = 0; i < n; i++) {
for (int j = JdkMath.max(i - 1, 0); j < n; j++) {
norm += JdkMath.abs(matrixT[i][j]);
}
}
// we can not handle a matrix with zero norm
if (Precision.equals(norm, 0.0, EPSILON)) {
throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
}
// Backsubstitute to find vectors of upper triangular form
double r = 0.0;
double s = 0.0;
double z = 0.0;
for (int idx = n - 1; idx >= 0; idx--) {
double p = realEigenvalues[idx];
double q = imagEigenvalues[idx];
if (Precision.equals(q, 0.0)) {
// Real vector
int l = idx;
matrixT[idx][idx] = 1.0;
for (int i = idx - 1; i >= 0; i--) {
double w = matrixT[i][i] - p;
r = 0.0;
for (int j = l; j <= idx; j++) {
r += matrixT[i][j] * matrixT[j][idx];
}
if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
z = w;
s = r;
} else {
l = i;
if (Precision.equals(imagEigenvalues[i], 0.0)) {
if (w != 0.0) {
matrixT[i][idx] = -r / w;
} else {
matrixT[i][idx] = -r / (Precision.EPSILON * norm);
}
} else {
// Solve real equations
double x = matrixT[i][i + 1];
double y = matrixT[i + 1][i];
q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
imagEigenvalues[i] * imagEigenvalues[i];
double t = (x * s - z * r) / q;
matrixT[i][idx] = t;
if (JdkMath.abs(x) > JdkMath.abs(z)) {
matrixT[i + 1][idx] = (-r - w * t) / x;
} else {
matrixT[i + 1][idx] = (-s - y * t) / z;
}
}
// Overflow control
double t = JdkMath.abs(matrixT[i][idx]);
if ((Precision.EPSILON * t) * t > 1) {
for (int j = i; j <= idx; j++) {
matrixT[j][idx] /= t;
}
}
}
}
} else if (q < 0.0) {
// Complex vector
int l = idx - 1;
// Last vector component imaginary so matrix is triangular
if (JdkMath.abs(matrixT[idx][idx - 1]) > JdkMath.abs(matrixT[idx - 1][idx])) {
matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
matrixT[idx - 1][idx] = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
} else {
final Complex result = cdiv(0.0, -matrixT[idx - 1][idx],
matrixT[idx - 1][idx - 1] - p, q);
matrixT[idx - 1][idx - 1] = result.getReal();
matrixT[idx - 1][idx] = result.getImaginary();
}
matrixT[idx][idx - 1] = 0.0;
matrixT[idx][idx] = 1.0;
for (int i = idx - 2; i >= 0; i--) {
double ra = 0.0;
double sa = 0.0;
for (int j = l; j <= idx; j++) {
ra += matrixT[i][j] * matrixT[j][idx - 1];
sa += matrixT[i][j] * matrixT[j][idx];
}
double w = matrixT[i][i] - p;
if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
z = w;
r = ra;
s = sa;
} else {
l = i;
if (Precision.equals(imagEigenvalues[i], 0.0)) {
final Complex c = cdiv(-ra, -sa, w, q);
matrixT[i][idx - 1] = c.getReal();
matrixT[i][idx] = c.getImaginary();
} else {
// Solve complex equations
double x = matrixT[i][i + 1];
double y = matrixT[i + 1][i];
double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
imagEigenvalues[i] * imagEigenvalues[i] - q * q;
final double vi = (realEigenvalues[i] - p) * 2.0 * q;
if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
vr = Precision.EPSILON * norm *
(JdkMath.abs(w) + JdkMath.abs(q) + JdkMath.abs(x) +
JdkMath.abs(y) + JdkMath.abs(z));
}
final Complex c = cdiv(x * r - z * ra + q * sa,
x * s - z * sa - q * ra, vr, vi);
matrixT[i][idx - 1] = c.getReal();
matrixT[i][idx] = c.getImaginary();
if (JdkMath.abs(x) > (JdkMath.abs(z) + JdkMath.abs(q))) {
matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] +
q * matrixT[i][idx]) / x;
matrixT[i + 1][idx] = (-sa - w * matrixT[i][idx] -
q * matrixT[i][idx - 1]) / x;
} else {
final Complex c2 = cdiv(-r - y * matrixT[i][idx - 1],
-s - y * matrixT[i][idx], z, q);
matrixT[i + 1][idx - 1] = c2.getReal();
matrixT[i + 1][idx] = c2.getImaginary();
}
}
// Overflow control
double t = JdkMath.max(JdkMath.abs(matrixT[i][idx - 1]),
JdkMath.abs(matrixT[i][idx]));
if ((Precision.EPSILON * t) * t > 1) {
for (int j = i; j <= idx; j++) {
matrixT[j][idx - 1] /= t;
matrixT[j][idx] /= t;
}
}
}
}
}
}
// Back transformation to get eigenvectors of original matrix
for (int j = n - 1; j >= 0; j--) {
for (int i = 0; i <= n - 1; i++) {
z = 0.0;
for (int k = 0; k <= JdkMath.min(j, n - 1); k++) {
z += matrixP[i][k] * matrixT[k][j];
}
matrixP[i][j] = z;
}
}
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] = matrixP[j][i];
}
eigenvectors[i] = new ArrayRealVector(tmp);
}
}