List of usage examples for org.apache.commons.math3.linear SchurTransformer getP
public RealMatrix getP()
From source file:org.pmad.gmm.MyEDC.java
/** * Find eigenvectors from a matrix transformed to Schur form. * * @param schur the schur transformation of the matrix * @throws MathArithmeticException if the Schur form has a norm of zero *///from w ww . ja va 2 s. c o m 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 = FastMath.max(i - 1, 0); j < n; j++) { norm += FastMath.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 (FastMath.abs(x) > FastMath.abs(z)) { matrixT[i + 1][idx] = (-r - w * t) / x; } else { matrixT[i + 1][idx] = (-s - y * t) / z; } } // Overflow control double t = FastMath.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 (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.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 * (FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) + FastMath.abs(y) + FastMath.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 (FastMath.abs(x) > (FastMath.abs(z) + FastMath.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 = FastMath.max(FastMath.abs(matrixT[i][idx - 1]), FastMath.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 <= FastMath.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); } }