A m × m matrix A can be written as the product of three matrices: A = P > + * × H × PT with P an unitary matrix and H a Hessenberg > + * matrix. Both P and H are m × m matrices.

> + *

Transformation to Hessenberg form is often not a goal by itself, but it is an > + * intermediate step in more general decomposition algorithms like > + * {@link EigenDecomposition eigen decomposition}. This class is therefore > + * intended for internal use by the library and is not public. As a consequence > + * of this explicitly limited scope, many methods directly returns references to > + * internal arrays, not copies.

> + *

This class is based on the method orthes in class EigenvalueDecomposition > + * from the JAMA library.

> + * > + * @see MathWorld > + * @see Householder Transformations > + * @version \$Id\$ > + * @since 3.1 > + */ > +class HessenbergTransformer { > + /** Householder vectors. */ > + private final double householderVectors[][]; > + /** Temporary storage vector. */ > + private final double ort[]; > + /** Cached value of P. */ > + private RealMatrix cachedP; > + /** Cached value of Pt. */ > + private RealMatrix cachedPt; > + /** Cached value of H. */ > + private RealMatrix cachedH; > + > + /** > + * Build the transformation to Hessenberg form of a general matrix. > + * > + * @param matrix matrix to transform. > + * @throws NonSquareMatrixException if the matrix is not square. > + */ > + public HessenbergTransformer(RealMatrix matrix) { > + if (!matrix.isSquare()) { > + throw new NonSquareMatrixException(matrix.getRowDimension(), > + matrix.getColumnDimension()); > + } > + > + final int m = matrix.getRowDimension(); > + householderVectors = matrix.getData(); > + ort = new double[m]; > + cachedP = null; > + cachedPt = null; > + cachedH = null; > + > + // transform matrix > + transform(); > + } > + > + /** > + * Returns the matrix P of the transform. > + *

P is an unitary matrix, i.e. its inverse is also its transpose.

> + * > + * @return the P matrix > + */ > + public RealMatrix getP() { > + if (cachedP == null) { > + final int n = householderVectors.length; > + final int high = n - 1; > + final double[][] pa = new double[n][n]; > + > + for (int i = 0; i < n; i++) { > + for (int j = 0; j < n; j++) { > + pa[i][j] = (i == j) ? 1 : 0; > + } > + } > + > + for (int m = high - 1; m >= 1; m--) { > + if (householderVectors[m][m - 1] != 0.0) { > + for (int i = m + 1; i <= high; i++) { > + ort[i] = householderVectors[i][m - 1]; > + } > + > + for (int j = m; j <= high; j++) { > + double g = 0.0; > + > + for (int i = m; i <= high; i++) { > + g += ort[i] * pa[i][j]; > + } > + > + // Double division avoids possible underflow > + g = (g / ort[m]) / householderVectors[m][m - 1]; > + > + for (int i = m; i <= high; i++) { > + pa[i][j] += g * ort[i]; > + } > + } > + } > + } > + > + cachedP = MatrixUtils.createRealMatrix(pa); > + } > + return cachedP; > + } > + > + /** > + * Returns the transpose of the matrix P of the transform. > + *

P is an unitary matrix, i.e. its inverse is also its transpose.

> + * > + * @return the transpose of the P matrix > + */ > + public RealMatrix getPT() { > + if (cachedPt == null) { > + cachedPt = getP().transpose(); > + } > + > + // return the cached matrix > + return cachedPt; > + } > + > + /** > + * Returns the Hessenberg matrix H of the transform. > + * > + * @return the H matrix > + */ > + public RealMatrix getH() { > + if (cachedH == null) { > + final int m = householderVectors.length; > + final double[][] h = new double[m][m]; > + for (int i = 0; i < m; ++i) { > + if (i > 0) { > + // copy the entry of the lower sub-diagonal > + h[i][i - 1] = householderVectors[i][i - 1]; > + } > + > + // copy upper triangular part of the matrix > + for (int j = i; j < m; ++j) { > + h[i][j] = householderVectors[i][j]; > + } > + } > + cachedH = MatrixUtils.createRealMatrix(h); > + } > + > + // return the cached matrix > + return cachedH; > + } > + > + /** > + * Get the Householder vectors of the transform. > + *

Note that since this class is only intended for internal use, it returns > + * directly a reference to its internal arrays, not a copy.

> + * > + * @return the main diagonal elements of the B matrix > + */ > + double[][] getHouseholderVectorsRef() { > + return householderVectors; > + } > + > + /** > + * Transform original matrix to Hessenberg form. > + *

Transformation is done using Householder transforms.