commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/SchurTransformer.java - commons-math - Git at Google

 /*
  * 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/LICENSE-2.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.math4.legacy.linear;

 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
 import org.apache.commons.math4.core.jdkmath.JdkMath;
 import org.apache.commons.numbers.core.Precision;

 /**
  * Class transforming a general real matrix to Schur form.
  * <p>A m &times; m matrix A can be written as the product of three matrices: A = P
  * &times; T &times; P<sup>T</sup> with P an orthogonal matrix and T an quasi-triangular
  * matrix. Both P and T are m &times; m matrices.</p>
  * <p>Transformation to Schur 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 hqr2 in class EigenvalueDecomposition
  * from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
  *
  * @see <a href="http://mathworld.wolfram.com/SchurDecomposition.html">Schur Decomposition - MathWorld</a>
  * @see <a href="http://en.wikipedia.org/wiki/Schur_decomposition">Schur Decomposition - Wikipedia</a>
  * @see <a href="http://en.wikipedia.org/wiki/Householder_transformation">Householder Transformations</a>
  * @since 3.1
  */
 class SchurTransformer {
     /** Maximum allowed iterations for convergence of the transformation. */
     private static final int MAX_ITERATIONS = 100;

     /** P matrix. */
     private final double[][] matrixP;
     /** T matrix. */
     private final double[][] matrixT;
     /** Cached value of P. */
     private RealMatrix cachedP;
     /** Cached value of T. */
     private RealMatrix cachedT;
     /** Cached value of PT. */
     private RealMatrix cachedPt;

     /** Epsilon criteria taken from JAMA code (originally was 2^-52). */
     private final double epsilon = Precision.EPSILON;

     /**
      * Build the transformation to Schur form of a general real matrix.
      *
      * @param matrix matrix to transform
      * @throws NonSquareMatrixException if the matrix is not square
      */
     SchurTransformer(final RealMatrix matrix) {
         if (!matrix.isSquare()) {
             throw new NonSquareMatrixException(matrix.getRowDimension(),
                                                matrix.getColumnDimension());
         }

         HessenbergTransformer transformer = new HessenbergTransformer(matrix);
         matrixT = transformer.getH().getData();
         matrixP = transformer.getP().getData();
         cachedT = null;
         cachedP = null;
         cachedPt = 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.</p>
      *
      * @return the P matrix
      */
     public RealMatrix getP() {
         if (cachedP == null) {
             cachedP = MatrixUtils.createRealMatrix(matrixP);
         }
         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.</p>
      *
      * @return the transpose of the P matrix
      */
     public RealMatrix getPT() {
         if (cachedPt == null) {
             cachedPt = getP().transpose();
         }

         // return the cached matrix
         return cachedPt;
     }

     /**
      * Returns the quasi-triangular Schur matrix T of the transform.
      *
      * @return the T matrix
      */
     public RealMatrix getT() {
         if (cachedT == null) {
             cachedT = MatrixUtils.createRealMatrix(matrixT);
         }

         // return the cached matrix
         return cachedT;
     }

     /**
      * Transform original matrix to Schur form.
      * @throws MaxCountExceededException if the transformation does not converge
      */
     private void transform() {
         final int n = matrixT.length;

         // compute matrix norm
         final double norm = getNorm();

         // shift information
         final ShiftInfo shift = new ShiftInfo();

         // Outer loop over eigenvalue index
         int iteration = 0;
         int iu = n - 1;
         while (iu >= 0) {

             // Look for single small sub-diagonal element
             final int il = findSmallSubDiagonalElement(iu, norm);

             // Check for convergence
             if (il == iu) {
                 // One root found
                 matrixT[iu][iu] += shift.exShift;
                 iu--;
                 iteration = 0;
             } else if (il == iu - 1) {
                 // Two roots found
                 double p = (matrixT[iu - 1][iu - 1] - matrixT[iu][iu]) / 2.0;
                 double q = p * p + matrixT[iu][iu - 1] * matrixT[iu - 1][iu];
                 matrixT[iu][iu] += shift.exShift;
                 matrixT[iu - 1][iu - 1] += shift.exShift;

                 if (q >= 0) {
                     double z = JdkMath.sqrt(JdkMath.abs(q));
                     if (p >= 0) {
                         z = p + z;
                     } else {
                         z = p - z;
                     }
                     final double x = matrixT[iu][iu - 1];
                     final double s = JdkMath.abs(x) + JdkMath.abs(z);
                     p = x / s;
                     q = z / s;
                     final double r = JdkMath.sqrt(p * p + q * q);
                     p /= r;
                     q /= r;

                     // Row modification
                     for (int j = iu - 1; j < n; j++) {
                         z = matrixT[iu - 1][j];
                         matrixT[iu - 1][j] = q * z + p * matrixT[iu][j];
                         matrixT[iu][j] = q * matrixT[iu][j] - p * z;
                     }

                     // Column modification
                     for (int i = 0; i <= iu; i++) {
                         z = matrixT[i][iu - 1];
                         matrixT[i][iu - 1] = q * z + p * matrixT[i][iu];
                         matrixT[i][iu] = q * matrixT[i][iu] - p * z;
                     }

                     // Accumulate transformations
                     for (int i = 0; i <= n - 1; i++) {
                         z = matrixP[i][iu - 1];
                         matrixP[i][iu - 1] = q * z + p * matrixP[i][iu];
                         matrixP[i][iu] = q * matrixP[i][iu] - p * z;
                     }
                 }
                 iu -= 2;
                 iteration = 0;
             } else {
                 // No convergence yet
                 computeShift(il, iu, iteration, shift);

                 // stop transformation after too many iterations
                 if (++iteration > MAX_ITERATIONS) {
                     throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
                                                         MAX_ITERATIONS);
                 }

                 // the initial houseHolder vector for the QR step
                 final double[] hVec = new double[3];

                 final int im = initQRStep(il, iu, shift, hVec);
                 performDoubleQRStep(il, im, iu, shift, hVec);
             }
         }
     }

     /**
      * Computes the L1 norm of the (quasi-)triangular matrix T.
      *
      * @return the L1 norm of matrix T
      */
     private double getNorm() {
         double norm = 0.0;
         for (int i = 0; i < matrixT.length; i++) {
             // as matrix T is (quasi-)triangular, also take the sub-diagonal element into account
             for (int j = JdkMath.max(i - 1, 0); j < matrixT.length; j++) {
                 norm += JdkMath.abs(matrixT[i][j]);
             }
         }
         return norm;
     }

     /**
      * Find the first small sub-diagonal element and returns its index.
      *
      * @param startIdx the starting index for the search
      * @param norm the L1 norm of the matrix
      * @return the index of the first small sub-diagonal element
      */
     private int findSmallSubDiagonalElement(final int startIdx, final double norm) {
         int l = startIdx;
         while (l > 0) {
             double s = JdkMath.abs(matrixT[l - 1][l - 1]) + JdkMath.abs(matrixT[l][l]);
             if (s == 0.0) {
                 s = norm;
             }
             if (JdkMath.abs(matrixT[l][l - 1]) < epsilon * s) {
                 break;
             }
             l--;
         }
         return l;
     }

     /**
      * Compute the shift for the current iteration.
      *
      * @param l the index of the small sub-diagonal element
      * @param idx the current eigenvalue index
      * @param iteration the current iteration
      * @param shift holder for shift information
      */
     private void computeShift(final int l, final int idx, final int iteration, final ShiftInfo shift) {
         // Form shift
         shift.x = matrixT[idx][idx];
         shift.y = shift.w = 0.0;
         if (l < idx) {
             shift.y = matrixT[idx - 1][idx - 1];
             shift.w = matrixT[idx][idx - 1] * matrixT[idx - 1][idx];
         }

         // Wilkinson's original ad hoc shift
         if (iteration == 10) {
             shift.exShift += shift.x;
             for (int i = 0; i <= idx; i++) {
                 matrixT[i][i] -= shift.x;
             }
             final double s = JdkMath.abs(matrixT[idx][idx - 1]) + JdkMath.abs(matrixT[idx - 1][idx - 2]);
             shift.x = 0.75 * s;
             shift.y = 0.75 * s;
             shift.w = -0.4375 * s * s;
         }

         // MATLAB's new ad hoc shift
         if (iteration == 30) {
             double s = (shift.y - shift.x) / 2.0;
             s = s * s + shift.w;
             if (s > 0.0) {
                 s = JdkMath.sqrt(s);
                 if (shift.y < shift.x) {
                     s = -s;
                 }
                 s = shift.x - shift.w / ((shift.y - shift.x) / 2.0 + s);
                 for (int i = 0; i <= idx; i++) {
                     matrixT[i][i] -= s;
                 }
                 shift.exShift += s;
                 shift.x = shift.y = shift.w = 0.964;
             }
         }
     }

     /**
      * Initialize the householder vectors for the QR step.
      *
      * @param il the index of the small sub-diagonal element
      * @param iu the current eigenvalue index
      * @param shift shift information holder
      * @param hVec the initial houseHolder vector
      * @return the start index for the QR step
      */
     private int initQRStep(int il, final int iu, final ShiftInfo shift, double[] hVec) {
         // Look for two consecutive small sub-diagonal elements
         int im = iu - 2;
         while (im >= il) {
             final double z = matrixT[im][im];
             final double r = shift.x - z;
             double s = shift.y - z;
             hVec[0] = (r * s - shift.w) / matrixT[im + 1][im] + matrixT[im][im + 1];
             hVec[1] = matrixT[im + 1][im + 1] - z - r - s;
             hVec[2] = matrixT[im + 2][im + 1];

             if (im == il) {
                 break;
             }

             final double lhs = JdkMath.abs(matrixT[im][im - 1]) * (JdkMath.abs(hVec[1]) + JdkMath.abs(hVec[2]));
             final double rhs = JdkMath.abs(hVec[0]) * (JdkMath.abs(matrixT[im - 1][im - 1]) +
                                                         JdkMath.abs(z) +
                                                         JdkMath.abs(matrixT[im + 1][im + 1]));

             if (lhs < epsilon * rhs) {
                 break;
             }
             im--;
         }

         return im;
     }

     /**
      * Perform a double QR step involving rows l:idx and columns m:n.
      *
      * @param il the index of the small sub-diagonal element
      * @param im the start index for the QR step
      * @param iu the current eigenvalue index
      * @param shift shift information holder
      * @param hVec the initial houseHolder vector
      */
     private void performDoubleQRStep(final int il, final int im, final int iu,
                                      final ShiftInfo shift, final double[] hVec) {

         final int n = matrixT.length;
         double p = hVec[0];
         double q = hVec[1];
         double r = hVec[2];

         for (int k = im; k <= iu - 1; k++) {
             boolean notlast = k != (iu - 1);
             if (k != im) {
                 p = matrixT[k][k - 1];
                 q = matrixT[k + 1][k - 1];
                 r = notlast ? matrixT[k + 2][k - 1] : 0.0;
                 shift.x = JdkMath.abs(p) + JdkMath.abs(q) + JdkMath.abs(r);
                 if (Precision.equals(shift.x, 0.0, epsilon)) {
                     continue;
                 }
                 p /= shift.x;
                 q /= shift.x;
                 r /= shift.x;
             }
             double s = JdkMath.sqrt(p * p + q * q + r * r);
             if (p < 0.0) {
                 s = -s;
             }
             if (s != 0.0) {
                 if (k != im) {
                     matrixT[k][k - 1] = -s * shift.x;
                 } else if (il != im) {
                     matrixT[k][k - 1] = -matrixT[k][k - 1];
                 }
                 p += s;
                 shift.x = p / s;
                 shift.y = q / s;
                 double z = r / s;
                 q /= p;
                 r /= p;

                 // Row modification
                 for (int j = k; j < n; j++) {
                     p = matrixT[k][j] + q * matrixT[k + 1][j];
                     if (notlast) {
                         p += r * matrixT[k + 2][j];
                         matrixT[k + 2][j] -= p * z;
                     }
                     matrixT[k][j] -= p * shift.x;
                     matrixT[k + 1][j] -= p * shift.y;
                 }

                 // Column modification
                 for (int i = 0; i <= JdkMath.min(iu, k + 3); i++) {
                     p = shift.x * matrixT[i][k] + shift.y * matrixT[i][k + 1];
                     if (notlast) {
                         p += z * matrixT[i][k + 2];
                         matrixT[i][k + 2] -= p * r;
                     }
                     matrixT[i][k] -= p;
                     matrixT[i][k + 1] -= p * q;
                 }

                 // Accumulate transformations
                 final int high = matrixT.length - 1;
                 for (int i = 0; i <= high; i++) {
                     p = shift.x * matrixP[i][k] + shift.y * matrixP[i][k + 1];
                     if (notlast) {
                         p += z * matrixP[i][k + 2];
                         matrixP[i][k + 2] -= p * r;
                     }
                     matrixP[i][k] -= p;
                     matrixP[i][k + 1] -= p * q;
                 }
             }  // (s != 0)
         }  // k loop

         // clean up pollution due to round-off errors
         for (int i = im + 2; i <= iu; i++) {
             matrixT[i][i-2] = 0.0;
             if (i > im + 2) {
                 matrixT[i][i-3] = 0.0;
             }
         }
     }

     /**
      * Internal data structure holding the current shift information.
      * Contains variable names as present in the original JAMA code.
      */
     private static class ShiftInfo {
         // CHECKSTYLE: stop all

         /** x shift info. */
         double x;
         /** y shift info. */
         double y;
         /** w shift info. */
         double w;
         /** Indicates an exceptional shift. */
         double exShift;

         // CHECKSTYLE: resume all
     }
 }
	/*
	* 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/LICENSE-2.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.math4.legacy.linear;

	import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
	import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
	import org.apache.commons.math4.core.jdkmath.JdkMath;
	import org.apache.commons.numbers.core.Precision;

	/**
	* Class transforming a general real matrix to Schur form.
	* <p>A m × m matrix A can be written as the product of three matrices: A = P
	* × T × P<sup>T</sup> with P an orthogonal matrix and T an quasi-triangular
	* matrix. Both P and T are m × m matrices.</p>
	* <p>Transformation to Schur 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 hqr2 in class EigenvalueDecomposition
	* from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
	*
	* @see <a href="http://mathworld.wolfram.com/SchurDecomposition.html">Schur Decomposition - MathWorld</a>
	* @see <a href="http://en.wikipedia.org/wiki/Schur_decomposition">Schur Decomposition - Wikipedia</a>
	* @see <a href="http://en.wikipedia.org/wiki/Householder_transformation">Householder Transformations</a>
	* @since 3.1
	*/
	class SchurTransformer {
	/** Maximum allowed iterations for convergence of the transformation. */
	private static final int MAX_ITERATIONS = 100;

	/** P matrix. */
	private final double[][] matrixP;
	/** T matrix. */
	private final double[][] matrixT;
	/** Cached value of P. */
	private RealMatrix cachedP;
	/** Cached value of T. */
	private RealMatrix cachedT;
	/** Cached value of PT. */
	private RealMatrix cachedPt;

	/** Epsilon criteria taken from JAMA code (originally was 2^-52). */
	private final double epsilon = Precision.EPSILON;

	/**
	* Build the transformation to Schur form of a general real matrix.
	*
	* @param matrix matrix to transform
	* @throws NonSquareMatrixException if the matrix is not square
	*/
	SchurTransformer(final RealMatrix matrix) {
	if (!matrix.isSquare()) {
	throw new NonSquareMatrixException(matrix.getRowDimension(),
	matrix.getColumnDimension());
	}

	HessenbergTransformer transformer = new HessenbergTransformer(matrix);
	matrixT = transformer.getH().getData();
	matrixP = transformer.getP().getData();
	cachedT = null;
	cachedP = null;
	cachedPt = 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.</p>
	*
	* @return the P matrix
	*/
	public RealMatrix getP() {
	if (cachedP == null) {
	cachedP = MatrixUtils.createRealMatrix(matrixP);
	}
	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.</p>
	*
	* @return the transpose of the P matrix
	*/
	public RealMatrix getPT() {
	if (cachedPt == null) {
	cachedPt = getP().transpose();
	}

	// return the cached matrix
	return cachedPt;
	}

	/**
	* Returns the quasi-triangular Schur matrix T of the transform.
	*
	* @return the T matrix
	*/
	public RealMatrix getT() {
	if (cachedT == null) {
	cachedT = MatrixUtils.createRealMatrix(matrixT);
	}

	// return the cached matrix
	return cachedT;
	}

	/**
	* Transform original matrix to Schur form.
	* @throws MaxCountExceededException if the transformation does not converge
	*/
	private void transform() {
	final int n = matrixT.length;

	// compute matrix norm
	final double norm = getNorm();

	// shift information
	final ShiftInfo shift = new ShiftInfo();

	// Outer loop over eigenvalue index
	int iteration = 0;
	int iu = n - 1;
	while (iu >= 0) {

	// Look for single small sub-diagonal element
	final int il = findSmallSubDiagonalElement(iu, norm);

	// Check for convergence
	if (il == iu) {
	// One root found
	matrixT[iu][iu] += shift.exShift;
	iu--;
	iteration = 0;
	} else if (il == iu - 1) {
	// Two roots found
	double p = (matrixT[iu - 1][iu - 1] - matrixT[iu][iu]) / 2.0;
	double q = p * p + matrixT[iu][iu - 1] * matrixT[iu - 1][iu];
	matrixT[iu][iu] += shift.exShift;
	matrixT[iu - 1][iu - 1] += shift.exShift;

	if (q >= 0) {
	double z = JdkMath.sqrt(JdkMath.abs(q));
	if (p >= 0) {
	z = p + z;
	} else {
	z = p - z;
	}
	final double x = matrixT[iu][iu - 1];
	final double s = JdkMath.abs(x) + JdkMath.abs(z);
	p = x / s;
	q = z / s;
	final double r = JdkMath.sqrt(p * p + q * q);
	p /= r;
	q /= r;

	// Row modification
	for (int j = iu - 1; j < n; j++) {
	z = matrixT[iu - 1][j];
	matrixT[iu - 1][j] = q * z + p * matrixT[iu][j];
	matrixT[iu][j] = q * matrixT[iu][j] - p * z;
	}

	// Column modification
	for (int i = 0; i <= iu; i++) {
	z = matrixT[i][iu - 1];
	matrixT[i][iu - 1] = q * z + p * matrixT[i][iu];
	matrixT[i][iu] = q * matrixT[i][iu] - p * z;
	}

	// Accumulate transformations
	for (int i = 0; i <= n - 1; i++) {
	z = matrixP[i][iu - 1];
	matrixP[i][iu - 1] = q * z + p * matrixP[i][iu];
	matrixP[i][iu] = q * matrixP[i][iu] - p * z;
	}
	}
	iu -= 2;
	iteration = 0;
	} else {
	// No convergence yet
	computeShift(il, iu, iteration, shift);

	// stop transformation after too many iterations
	if (++iteration > MAX_ITERATIONS) {
	throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
	MAX_ITERATIONS);
	}

	// the initial houseHolder vector for the QR step
	final double[] hVec = new double[3];

	final int im = initQRStep(il, iu, shift, hVec);
	performDoubleQRStep(il, im, iu, shift, hVec);
	}
	}
	}

	/**
	* Computes the L1 norm of the (quasi-)triangular matrix T.
	*
	* @return the L1 norm of matrix T
	*/
	private double getNorm() {
	double norm = 0.0;
	for (int i = 0; i < matrixT.length; i++) {
	// as matrix T is (quasi-)triangular, also take the sub-diagonal element into account
	for (int j = JdkMath.max(i - 1, 0); j < matrixT.length; j++) {
	norm += JdkMath.abs(matrixT[i][j]);
	}
	}
	return norm;
	}

	/**
	* Find the first small sub-diagonal element and returns its index.
	*
	* @param startIdx the starting index for the search
	* @param norm the L1 norm of the matrix
	* @return the index of the first small sub-diagonal element
	*/
	private int findSmallSubDiagonalElement(final int startIdx, final double norm) {
	int l = startIdx;
	while (l > 0) {
	double s = JdkMath.abs(matrixT[l - 1][l - 1]) + JdkMath.abs(matrixT[l][l]);
	if (s == 0.0) {
	s = norm;
	}
	if (JdkMath.abs(matrixT[l][l - 1]) < epsilon * s) {
	break;
	}
	l--;
	}
	return l;
	}

	/**
	* Compute the shift for the current iteration.
	*
	* @param l the index of the small sub-diagonal element
	* @param idx the current eigenvalue index
	* @param iteration the current iteration
	* @param shift holder for shift information
	*/
	private void computeShift(final int l, final int idx, final int iteration, final ShiftInfo shift) {
	// Form shift
	shift.x = matrixT[idx][idx];
	shift.y = shift.w = 0.0;
	if (l < idx) {
	shift.y = matrixT[idx - 1][idx - 1];
	shift.w = matrixT[idx][idx - 1] * matrixT[idx - 1][idx];
	}

	// Wilkinson's original ad hoc shift
	if (iteration == 10) {
	shift.exShift += shift.x;
	for (int i = 0; i <= idx; i++) {
	matrixT[i][i] -= shift.x;
	}
	final double s = JdkMath.abs(matrixT[idx][idx - 1]) + JdkMath.abs(matrixT[idx - 1][idx - 2]);
	shift.x = 0.75 * s;
	shift.y = 0.75 * s;
	shift.w = -0.4375 * s * s;
	}

	// MATLAB's new ad hoc shift
	if (iteration == 30) {
	double s = (shift.y - shift.x) / 2.0;
	s = s * s + shift.w;
	if (s > 0.0) {
	s = JdkMath.sqrt(s);
	if (shift.y < shift.x) {
	s = -s;
	}
	s = shift.x - shift.w / ((shift.y - shift.x) / 2.0 + s);
	for (int i = 0; i <= idx; i++) {
	matrixT[i][i] -= s;
	}
	shift.exShift += s;
	shift.x = shift.y = shift.w = 0.964;
	}
	}
	}

	/**
	* Initialize the householder vectors for the QR step.
	*
	* @param il the index of the small sub-diagonal element
	* @param iu the current eigenvalue index
	* @param shift shift information holder
	* @param hVec the initial houseHolder vector
	* @return the start index for the QR step
	*/
	private int initQRStep(int il, final int iu, final ShiftInfo shift, double[] hVec) {
	// Look for two consecutive small sub-diagonal elements
	int im = iu - 2;
	while (im >= il) {
	final double z = matrixT[im][im];
	final double r = shift.x - z;
	double s = shift.y - z;
	hVec[0] = (r * s - shift.w) / matrixT[im + 1][im] + matrixT[im][im + 1];
	hVec[1] = matrixT[im + 1][im + 1] - z - r - s;
	hVec[2] = matrixT[im + 2][im + 1];

	if (im == il) {
	break;
	}

	final double lhs = JdkMath.abs(matrixT[im][im - 1]) * (JdkMath.abs(hVec[1]) + JdkMath.abs(hVec[2]));
	final double rhs = JdkMath.abs(hVec[0]) * (JdkMath.abs(matrixT[im - 1][im - 1]) +
	JdkMath.abs(z) +
	JdkMath.abs(matrixT[im + 1][im + 1]));

	if (lhs < epsilon * rhs) {
	break;
	}
	im--;
	}

	return im;
	}

	/**
	* Perform a double QR step involving rows l:idx and columns m:n.
	*
	* @param il the index of the small sub-diagonal element
	* @param im the start index for the QR step
	* @param iu the current eigenvalue index
	* @param shift shift information holder
	* @param hVec the initial houseHolder vector
	*/
	private void performDoubleQRStep(final int il, final int im, final int iu,
	final ShiftInfo shift, final double[] hVec) {

	final int n = matrixT.length;
	double p = hVec[0];
	double q = hVec[1];
	double r = hVec[2];

	for (int k = im; k <= iu - 1; k++) {
	boolean notlast = k != (iu - 1);
	if (k != im) {
	p = matrixT[k][k - 1];
	q = matrixT[k + 1][k - 1];
	r = notlast ? matrixT[k + 2][k - 1] : 0.0;
	shift.x = JdkMath.abs(p) + JdkMath.abs(q) + JdkMath.abs(r);
	if (Precision.equals(shift.x, 0.0, epsilon)) {
	continue;
	}
	p /= shift.x;
	q /= shift.x;
	r /= shift.x;
	}
	double s = JdkMath.sqrt(p * p + q * q + r * r);
	if (p < 0.0) {
	s = -s;
	}
	if (s != 0.0) {
	if (k != im) {
	matrixT[k][k - 1] = -s * shift.x;
	} else if (il != im) {
	matrixT[k][k - 1] = -matrixT[k][k - 1];
	}
	p += s;
	shift.x = p / s;
	shift.y = q / s;
	double z = r / s;
	q /= p;
	r /= p;

	// Row modification
	for (int j = k; j < n; j++) {
	p = matrixT[k][j] + q * matrixT[k + 1][j];
	if (notlast) {
	p += r * matrixT[k + 2][j];
	matrixT[k + 2][j] -= p * z;
	}
	matrixT[k][j] -= p * shift.x;
	matrixT[k + 1][j] -= p * shift.y;
	}

	// Column modification
	for (int i = 0; i <= JdkMath.min(iu, k + 3); i++) {
	p = shift.x * matrixT[i][k] + shift.y * matrixT[i][k + 1];
	if (notlast) {
	p += z * matrixT[i][k + 2];
	matrixT[i][k + 2] -= p * r;
	}
	matrixT[i][k] -= p;
	matrixT[i][k + 1] -= p * q;
	}

	// Accumulate transformations
	final int high = matrixT.length - 1;
	for (int i = 0; i <= high; i++) {
	p = shift.x * matrixP[i][k] + shift.y * matrixP[i][k + 1];
	if (notlast) {
	p += z * matrixP[i][k + 2];
	matrixP[i][k + 2] -= p * r;
	}
	matrixP[i][k] -= p;
	matrixP[i][k + 1] -= p * q;
	}
	} // (s != 0)
	} // k loop

	// clean up pollution due to round-off errors
	for (int i = im + 2; i <= iu; i++) {
	matrixT[i][i-2] = 0.0;
	if (i > im + 2) {
	matrixT[i][i-3] = 0.0;
	}
	}
	}

	/**
	* Internal data structure holding the current shift information.
	* Contains variable names as present in the original JAMA code.
	*/
	private static class ShiftInfo {
	// CHECKSTYLE: stop all

	/** x shift info. */
	double x;
	/** y shift info. */
	double y;
	/** w shift info. */
	double w;
	/** Indicates an exceptional shift. */
	double exShift;

	// CHECKSTYLE: resume all
	}
	}