

Java example source code file (HessenbergTransformer.java)
The HessenbergTransformer.java Java example source code/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math3.linear; import org.apache.commons.math3.util.FastMath; import org.apache.commons.math3.util.Precision; /** * Class transforming a general real matrix to Hessenberg form. * <p>A m × m matrix A can be written as the product of three matrices: A = P * × H × P<sup>T with P an orthogonal matrix and H a Hessenberg * matrix. Both P and H are m × m matrices.</p> * <p>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.</p> * <p>This class is based on the method orthes in class EigenvalueDecomposition * from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA library. * * @see <a href="http://mathworld.wolfram.com/HessenbergDecomposition.html">MathWorld * @see <a href="http://en.wikipedia.org/wiki/Householder_transformation">Householder Transformations * @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 */ HessenbergTransformer(final 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>P is an orthogonal 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>P is an orthogonal 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 subdiagonal 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. * <p>Note that since this class is only intended for internal use, it returns * directly a reference to its internal arrays, not a copy.</p> * * @return the main diagonal elements of the B matrix */ double[][] getHouseholderVectorsRef() { return householderVectors; } /** * Transform original matrix to Hessenberg form. * <p>Transformation is done using Householder transforms. */ private void transform() { final int n = householderVectors.length; final int high = n  1; for (int m = 1; m <= high  1; m++) { // Scale column. double scale = 0; for (int i = m; i <= high; i++) { scale += FastMath.abs(householderVectors[i][m  1]); } if (!Precision.equals(scale, 0)) { // Compute Householder transformation. double h = 0; for (int i = high; i >= m; i) { ort[i] = householderVectors[i][m  1] / scale; h += ort[i] * ort[i]; } final double g = (ort[m] > 0) ? FastMath.sqrt(h) : FastMath.sqrt(h); h = ort[m] * g; ort[m] = g; // Apply Householder similarity transformation // H = (I  u*u' / h) * H * (I  u*u' / h) for (int j = m; j < n; j++) { double f = 0; for (int i = high; i >= m; i) { f += ort[i] * householderVectors[i][j]; } f /= h; for (int i = m; i <= high; i++) { householderVectors[i][j] = f * ort[i]; } } for (int i = 0; i <= high; i++) { double f = 0; for (int j = high; j >= m; j) { f += ort[j] * householderVectors[i][j]; } f /= h; for (int j = m; j <= high; j++) { householderVectors[i][j] = f * ort[j]; } } ort[m] = scale * ort[m]; householderVectors[m][m  1] = scale * g; } } } } Other Java examples (source code examples)Here is a short list of links related to this Java HessenbergTransformer.java source code file: 
Copyright 19982019 Alvin Alexander, alvinalexander.com
All Rights Reserved.
A percentage of advertising revenue from
pages under the /java/jwarehouse
URI on this website is
paid back to open source projects.