/*
 * 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.math.linear;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.util.MathUtils;

/**
 * Calculates the eigen decomposition of a <strong>symmetric</strong> matrix.
 * <p>The eigen decomposition of matrix A is a set of two matrices:
 * V and D such that A = V D V<sup>T</sup>. A, V and D are all m &times; m
 * matrices.</p>
 * <p>As of 2.0, this class supports only <strong>symmetric</strong> matrices,
 * and hence computes only real realEigenvalues. This implies the D matrix returned by
 * {@link #getD()} is always diagonal and the imaginary values returned {@link
 * #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always null.</p>
 * <p>When called with a {@link RealMatrix} argument, this implementation only uses
 * the upper part of the matrix, the part below the diagonal is not accessed at all.</p>
 * <p>Eigenvalues are computed as soon as the matrix is decomposed, but eigenvectors
 * are computed only when required, i.e. only when one of the {@link #getEigenvector(int)},
 * {@link #getV()}, {@link #getVT()}, {@link #getSolver()} methods is called.</p>
 * <p>This implementation is based on Inderjit Singh Dhillon thesis
 * <a href="http://www.cs.utexas.edu/users/inderjit/public_papers/thesis.pdf">A
 * New O(n<sup>2</sup>) Algorithm for the Symmetric Tridiagonal Eigenvalue/Eigenvector
 * Problem</a>, on Beresford N. Parlett and Osni A. Marques paper <a
 * href="http://www.netlib.org/lapack/lawnspdf/lawn155.pdf">An Implementation of the
 * dqds Algorithm (Positive Case)</a> and on the corresponding LAPACK routines (DLARRE,
 * DLASQ2, DLAZQ3, DLAZQ4, DLASQ5 and DLASQ6).</p>
 * @author Beresford Parlett, University of California, Berkeley, USA (fortran version)
 * @author Jim Demmel, University of California, Berkeley, USA (fortran version)
 * @author Inderjit Dhillon, University of Texas, Austin, USA(fortran version)
 * @author Osni Marques, LBNL/NERSC, USA (fortran version)
 * @author Christof Voemel, University of California, Berkeley, USA(fortran version)
 * @version $Revision$ $Date$
 * @since 2.0
 */
public class EigenDecompositionImpl implements EigenDecomposition {

    /** Tolerance. */
    private static final double TOLERANCE = 100 * MathUtils.EPSILON;

    /** Squared tolerance. */
    private static final double TOLERANCE_2 = TOLERANCE * TOLERANCE;

    /** Split tolerance. */
    private double splitTolerance;

    /** Main diagonal of the tridiagonal matrix. */
    private double[] main;

    /** Secondary diagonal of the tridiagonal matrix. */
    private double[] secondary;

    /** Squared secondary diagonal of the tridiagonal matrix. */
    private double[] squaredSecondary;

    /** Transformer to tridiagonal (may be null if matrix is already tridiagonal). */
    private TriDiagonalTransformer transformer;

    /** Lower bound of spectra. */
    private double lowerSpectra;

    /** Upper bound of spectra. */
    private double upperSpectra;

    /** Minimum pivot in the Sturm sequence. */
    private double minPivot;

    /** Current shift. */
    private double sigma;

    /** Low part of the current shift. */
    private double sigmaLow;

    /** Shift increment to apply. */
    private double tau;

    /** Work array for all decomposition algorithms. */
    private double[] work;

    /** Shift within qd array for ping-pong implementation. */
    private int pingPong;

    /** Max value of diagonal elements in current segment. */
    private double qMax;

    /** Min value of off-diagonal elements in current segment. */
    private double eMin;

    /** Type of the last dqds shift. */
    private int    tType;

    /** Minimal value on current state of the diagonal. */
    private double dMin;

    /** Minimal value on current state of the diagonal, excluding last element. */
    private double dMin1;

    /** Minimal value on current state of the diagonal, excluding last two elements. */
    private double dMin2;

    /** Last value on current state of the diagonal. */
    private double dN;

    /** Last but one value on current state of the diagonal. */
    private double dN1;

    /** Last but two on current state of the diagonal. */
    private double dN2;

    /** Shift ratio with respect to dMin used when tType == 6. */
    private double g;

    /** Real part of the realEigenvalues. */
    private double[] realEigenvalues;

    /** Imaginary part of the realEigenvalues. */
    private double[] imagEigenvalues;

    /** Eigenvectors. */
    private ArrayRealVector[] eigenvectors;

    /** Cached value of V. */
    private RealMatrix cachedV;

    /** Cached value of D. */
    private RealMatrix cachedD;

    /** Cached value of Vt. */
    private RealMatrix cachedVt;

    /**
     * Calculates the eigen decomposition of the given symmetric matrix. 
     * @param matrix The <strong>symmetric</strong> matrix to decompose.
     * @param splitTolerance tolerance on the off-diagonal elements relative to the
     * geometric mean to split the tridiagonal matrix (a suggested value is
     * {@link MathUtils#SAFE_MIN})
     * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
     * if algorithm fails to converge
     */
    public EigenDecompositionImpl(final RealMatrix matrix,
                                  final double splitTolerance)
        throws InvalidMatrixException {
        if (isSymmetric(matrix)) {
            this.splitTolerance = splitTolerance;
            transformToTridiagonal(matrix);
            decompose();
        } else {
            // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are NOT supported
            // see issue https://issues.apache.org/jira/browse/MATH-235
            throw new InvalidMatrixException("eigen decomposition of assymetric matrices not supported yet");
        }
    }

    /**
     * Calculates the eigen decomposition of the given tridiagonal symmetric matrix. 
     * @param main the main diagonal of the matrix (will be copied)
     * @param secondary the secondary diagonal of the matrix (will be copied)
     * @param splitTolerance tolerance on the off-diagonal elements relative to the
     * geometric mean to split the tridiagonal matrix (a suggested value is
     * {@link MathUtils#SAFE_MIN})
     * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
     * if algorithm fails to converge
     */
    public EigenDecompositionImpl(final double[] main, double[] secondary,
            final double splitTolerance)
        throws InvalidMatrixException {

        this.main      = main.clone();
        this.secondary = secondary.clone();
        transformer    = null;

        // pre-compute some elements
        squaredSecondary = new double[secondary.length];
        for (int i = 0; i < squaredSecondary.length; ++i) {
            final double s = secondary[i];
            squaredSecondary[i] = s * s;
        }

        this.splitTolerance = splitTolerance;
        decompose();

    }

    /**
     * Check if a matrix is symmetric.
     * @param matrix matrix to check
     * @return true if matrix is symmetric
     */
    private boolean isSymmetric(final RealMatrix matrix) {
        final int rows    = matrix.getRowDimension();
        final int columns = matrix.getColumnDimension();
        final double eps  = 10 * rows * columns * MathUtils.EPSILON;
        for (int i = 0; i < rows; ++i) {
            for (int j = i + 1; j < columns; ++j) {
                final double mij = matrix.getEntry(i, j);
                final double mji = matrix.getEntry(j, i);
                if (Math.abs(mij - mji) > (Math.max(Math.abs(mij), Math.abs(mji)) * eps)) {
                    return false;
                }
            }
        }
        return true;
    }

    /**
     * Decompose a tridiagonal symmetric matrix. 
     * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
     * if algorithm fails to converge
     */
    private void decompose() {

        cachedV  = null;
        cachedD  = null;
        cachedVt = null;
        work     = new double[6 * main.length];

        // compute the Gershgorin circles
        computeGershgorinCircles();

        // find all the realEigenvalues
        findEigenvalues();

        // we will search for eigenvectors only if required
        eigenvectors = null;

    }

    /** {@inheritDoc} */
    public RealMatrix getV()
        throws InvalidMatrixException {

        if (cachedV == null) {

            if (eigenvectors == null) {
                findEigenVectors();
            }

            final int m = eigenvectors.length;
            cachedV = MatrixUtils.createRealMatrix(m, m);
            for (int k = 0; k < m; ++k) {
                cachedV.setColumnVector(k, eigenvectors[k]);
            }

        }

        // return the cached matrix
        return cachedV;

    }

    /** {@inheritDoc} */
    public RealMatrix getD()
        throws InvalidMatrixException {
        if (cachedD == null) {
            // cache the matrix for subsequent calls
            cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
        }
        return cachedD;
    }

    /** {@inheritDoc} */
    public RealMatrix getVT()
        throws InvalidMatrixException {

        if (cachedVt == null) {

            if (eigenvectors == null) {
                findEigenVectors();
            }

            final int m = eigenvectors.length;
            cachedVt = MatrixUtils.createRealMatrix(m, m);
            for (int k = 0; k < m; ++k) {
                cachedVt.setRowVector(k, eigenvectors[k]);
            }

        }

        // return the cached matrix
        return cachedVt;

    }

    /** {@inheritDoc} */
    public double[] getRealEigenvalues()
        throws InvalidMatrixException {
        return realEigenvalues.clone();
    }

    /** {@inheritDoc} */
    public double getRealEigenvalue(final int i)
        throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
        return realEigenvalues[i];
    }

    /** {@inheritDoc} */
    public double[] getImagEigenvalues()
        throws InvalidMatrixException {
        return imagEigenvalues.clone();
    }

    /** {@inheritDoc} */
    public double getImagEigenvalue(final int i)
        throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
        return imagEigenvalues[i];
    }

    /** {@inheritDoc} */
    public RealVector getEigenvector(final int i)
        throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
        if (eigenvectors == null) {
            findEigenVectors();
        }
        return eigenvectors[i].copy();
    }

    /**
     * Return the determinant of the matrix
     * @return determinant of the matrix
     */
    public double getDeterminant() {
        double determinant = 1;
        for (double lambda : realEigenvalues) {
            determinant *= lambda;
        }
        return determinant;
    }

    /** {@inheritDoc} */
    public DecompositionSolver getSolver() {
        if (eigenvectors == null) {
            findEigenVectors();
        }
        return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
    }

    /** Specialized solver. */
    private static class Solver implements DecompositionSolver {
    
        /** Real part of the realEigenvalues. */
        private double[] realEigenvalues;

        /** Imaginary part of the realEigenvalues. */
        private double[] imagEigenvalues;

        /** Eigenvectors. */
        private final ArrayRealVector[] eigenvectors;

        /**
         * Build a solver from decomposed matrix.
         * @param realEigenvalues real parts of the eigenvalues
         * @param imagEigenvalues imaginary parts of the eigenvalues
         * @param eigenvectors eigenvectors
         */
        private Solver(final double[] realEigenvalues, final double[] imagEigenvalues,
                       final ArrayRealVector[] eigenvectors) {
            this.realEigenvalues = realEigenvalues;
            this.imagEigenvalues = imagEigenvalues;
            this.eigenvectors    = eigenvectors; 
        }

        /** Solve the linear equation A &times; X = B for symmetric matrices A.
         * <p>This method only find exact linear solutions, i.e. solutions for
         * which ||A &times; X - B|| is exactly 0.</p>
         * @param b right-hand side of the equation A &times; X = B
         * @return a vector X that minimizes the two norm of A &times; X - B
         * @exception IllegalArgumentException if matrices dimensions don't match
         * @exception InvalidMatrixException if decomposed matrix is singular
         */
        public double[] solve(final double[] b)
            throws IllegalArgumentException, InvalidMatrixException {

            if (!isNonSingular()) {
                throw new SingularMatrixException();
            }

            final int m = realEigenvalues.length;
            if (b.length != m) {
                throw MathRuntimeException.createIllegalArgumentException(
                        "vector length mismatch: got {0} but expected {1}",
                        b.length, m);
            }

            final double[] bp = new double[m];
            for (int i = 0; i < m; ++i) {
                final ArrayRealVector v = eigenvectors[i];
                final double[] vData = v.getDataRef();
                final double s = v.dotProduct(b) / realEigenvalues[i];
                for (int j = 0; j < m; ++j) {
                    bp[j] += s * vData[j];
                }
            }

            return bp;

        }

        /** Solve the linear equation A &times; X = B for symmetric matrices A.
         * <p>This method only find exact linear solutions, i.e. solutions for
         * which ||A &times; X - B|| is exactly 0.</p>
         * @param b right-hand side of the equation A &times; X = B
         * @return a vector X that minimizes the two norm of A &times; X - B
         * @exception IllegalArgumentException if matrices dimensions don't match
         * @exception InvalidMatrixException if decomposed matrix is singular
         */
        public RealVector solve(final RealVector b)
            throws IllegalArgumentException, InvalidMatrixException {

            if (!isNonSingular()) {
                throw new SingularMatrixException();
            }

            final int m = realEigenvalues.length;
            if (b.getDimension() != m) {
                throw MathRuntimeException.createIllegalArgumentException(
                        "vector length mismatch: got {0} but expected {1}",
                        b.getDimension(), m);
            }

            final double[] bp = new double[m];
            for (int i = 0; i < m; ++i) {
                final ArrayRealVector v = eigenvectors[i];
                final double[] vData = v.getDataRef();
                final double s = v.dotProduct(b) / realEigenvalues[i];
                for (int j = 0; j < m; ++j) {
                    bp[j] += s * vData[j];
                }
            }

            return new ArrayRealVector(bp, false);

        }

        /** Solve the linear equation A &times; X = B for symmetric matrices A.
         * <p>This method only find exact linear solutions, i.e. solutions for
         * which ||A &times; X - B|| is exactly 0.</p>
         * @param b right-hand side of the equation A &times; X = B
         * @return a matrix X that minimizes the two norm of A &times; X - B
         * @exception IllegalArgumentException if matrices dimensions don't match
         * @exception InvalidMatrixException if decomposed matrix is singular
         */
        public RealMatrix solve(final RealMatrix b)
            throws IllegalArgumentException, InvalidMatrixException {

            if (!isNonSingular()) {
                throw new SingularMatrixException();
            }

            final int m = realEigenvalues.length;
            if (b.getRowDimension() != m) {
                throw MathRuntimeException.createIllegalArgumentException(
                        "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
                        b.getRowDimension(), b.getColumnDimension(), m, "n");
            }

            final int nColB = b.getColumnDimension();
            final double[][] bp = new double[m][nColB];
            for (int k = 0; k < nColB; ++k) {
                for (int i = 0; i < m; ++i) {
                    final ArrayRealVector v = eigenvectors[i];
                    final double[] vData = v.getDataRef();
                    double s = 0;
                    for (int j = 0; j < m; ++j) {
                        s += v.getEntry(j) * b.getEntry(j, k);
                    }
                    s /= realEigenvalues[i];
                    for (int j = 0; j < m; ++j) {
                        bp[j][k] += s * vData[j];
                    }
                }
            }

            return MatrixUtils.createRealMatrix(bp);

        }

        /**
         * Check if the decomposed matrix is non-singular.
         * @return true if the decomposed matrix is non-singular
         */
        public boolean isNonSingular() {
            for (int i = 0; i < realEigenvalues.length; ++i) {
                if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) {
                    return false;
                }
            }
            return true;
        }

        /** Get the inverse of the decomposed matrix.
         * @return inverse matrix
         * @throws InvalidMatrixException if decomposed matrix is singular
         */
        public RealMatrix getInverse()
            throws InvalidMatrixException {

            if (!isNonSingular()) {
                throw new SingularMatrixException();
            }

            final int m = realEigenvalues.length;
            final double[][] invData = new double[m][m];

            for (int i = 0; i < m; ++i) {
                final double[] invI = invData[i];
                for (int j = 0; j < m; ++j) {
                    double invIJ = 0;
                    for (int k = 0; k < m; ++k) {
                        final double[] vK = eigenvectors[k].getDataRef();
                        invIJ += vK[i] * vK[j] / realEigenvalues[k];
                    }
                    invI[j] = invIJ;
                }
            }
            return MatrixUtils.createRealMatrix(invData);

        }

    }

    /**
     * Transform matrix to tridiagonal.
     * @param matrix matrix to transform
     */
    private void transformToTridiagonal(final RealMatrix matrix) {

        // transform the matrix to tridiagonal
        transformer = new TriDiagonalTransformer(matrix);
        main      = transformer.getMainDiagonalRef();
        secondary = transformer.getSecondaryDiagonalRef();

        // pre-compute some elements
        squaredSecondary = new double[secondary.length];
        for (int i = 0; i < squaredSecondary.length; ++i) {
            final double s = secondary[i];
            squaredSecondary[i] = s * s;
        }

    }

    /**
     * Compute the Gershgorin circles for all rows.
     */
    private void computeGershgorinCircles() {

        final int m     = main.length;
        final int lowerStart = 4 * m;
        final int upperStart = 5 * m;
        lowerSpectra = Double.POSITIVE_INFINITY;
        upperSpectra = Double.NEGATIVE_INFINITY;
        double eMax = 0;

        double eCurrent = 0;
        for (int i = 0; i < m - 1; ++i) {

            final double dCurrent = main[i];
            final double ePrevious = eCurrent;
            eCurrent = Math.abs(secondary[i]);
            eMax = Math.max(eMax, eCurrent);
            final double radius = ePrevious + eCurrent;

            final double lower = dCurrent - radius;
            work[lowerStart + i] = lower;
            lowerSpectra = Math.min(lowerSpectra, lower);

            final double upper = dCurrent + radius;
            work[upperStart + i] = upper;
            upperSpectra = Math.max(upperSpectra, upper);
            
        }

        final double dCurrent = main[m - 1];
        work[lowerStart + m - 1] = dCurrent - eCurrent;
        work[upperStart + m - 1] = dCurrent + eCurrent;
        minPivot = MathUtils.SAFE_MIN * Math.max(1.0, eMax * eMax);

    }

    /**
     * Find the realEigenvalues.
     * @exception InvalidMatrixException if a block cannot be diagonalized
     */
    private void findEigenvalues()
        throws InvalidMatrixException {

        // compute splitting points
        List<Integer> splitIndices = computeSplits();

        // find realEigenvalues in each block
        realEigenvalues = new double[main.length];
        imagEigenvalues = new double[main.length];
        int begin = 0;
        for (final int end : splitIndices) {
            final int n = end - begin;
            switch (n) {

            case 1:
                // apply dedicated method for dimension 1
                process1RowBlock(begin);
                break;

            case 2:
                // apply dedicated method for dimension 2
                process2RowsBlock(begin);
                break;

            case 3:
                // apply dedicated method for dimension 3
                process3RowsBlock(begin);
                break;

            default:

                // choose an initial shift for LDL<sup>T</sup> decomposition
                final double[] range       = eigenvaluesRange(begin, n);
                final double oneFourth     = 0.25 * (3 * range[0] + range[1]);
                final int oneFourthCount   = countEigenValues(oneFourth, begin, n);
                final double threeFourth   = 0.25 * (range[0] + 3 * range[1]);
                final int threeFourthCount = countEigenValues(threeFourth, begin, n);
                final boolean chooseLeft   = (oneFourthCount - 1) >= (n - threeFourthCount);
                final double lambda        = chooseLeft ? range[0] : range[1];

                tau = (range[1] - range[0]) * MathUtils.EPSILON * n + 2 * minPivot;

                // decompose T&lambda;I as LDL<sup>T</sup>
                ldlTDecomposition(lambda, begin, n);

                // apply general dqd/dqds method
                processGeneralBlock(n);

                // extract realEigenvalues
                if (chooseLeft) {
                    for (int i = 0; i < n; ++i) {
                        realEigenvalues[begin + i] = lambda + work[4 * i];
                    }
                } else {
                    for (int i = 0; i < n; ++i) {
                        realEigenvalues[begin + i] = lambda - work[4 * i];
                    }                    
                }

            }
            begin = end;
        }

        // sort the realEigenvalues in decreasing order
        Arrays.sort(realEigenvalues);
        for (int i = 0, j = realEigenvalues.length - 1; i < j; ++i, --j) {
            final double tmp = realEigenvalues[i];
            realEigenvalues[i] = realEigenvalues[j];
            realEigenvalues[j] = tmp;
        }

    }

    /**
     * Compute splitting points.
     * @return list of indices after matrix can be split
     */
    private List<Integer> computeSplits() {

        final List<Integer> list = new ArrayList<Integer>();

        // splitting preserving relative accuracy
        double absDCurrent = Math.abs(main[0]);
        for (int i = 0; i < secondary.length; ++i) {
            final double absDPrevious = absDCurrent;
            absDCurrent = Math.abs(main[i + 1]);
            final double max = splitTolerance * Math.sqrt(absDPrevious * absDCurrent);
            if (Math.abs(secondary[i]) <= max) {
                list.add(i + 1);
                secondary[i] = 0;
                squaredSecondary[i] = 0;
            }
        }

        list.add(secondary.length + 1);
        return list;

    }

    /**
     * Find eigenvalue in a block with 1 row.
     * <p>In low dimensions, we simply solve the characteristic polynomial.</p>
     * @param index index of the first row of the block
     */
    private void process1RowBlock(final int index) {
        realEigenvalues[index] = main[index];
    }

    /**
     * Find realEigenvalues in a block with 2 rows.
     * <p>In low dimensions, we simply solve the characteristic polynomial.</p>
     * @param index index of the first row of the block
     * @exception InvalidMatrixException if characteristic polynomial cannot be solved
     */
    private void process2RowsBlock(final int index)
        throws InvalidMatrixException {

        // the characteristic polynomial is
        // X^2 - (q0 + q1) X + q0 q1 - e1^2
        final double q0   = main[index];
        final double q1   = main[index + 1];
        final double e12  = squaredSecondary[index];

        final double s     = q0 + q1;
        final double p     = q0 * q1 - e12;
        final double delta = s * s - 4 * p;
        if (delta < 0) {
            throw new InvalidMatrixException("cannot solve degree {0} equation", 2);
        }

        final double largestRoot = 0.5 * (s + Math.sqrt(delta));
        realEigenvalues[index]     = largestRoot;
        realEigenvalues[index + 1] = p / largestRoot;

    }

    /**
     * Find realEigenvalues in a block with 3 rows.
     * <p>In low dimensions, we simply solve the characteristic polynomial.</p>
     * @param index index of the first row of the block
     * @exception InvalidMatrixException if diagonal elements are not positive
     */
    private void process3RowsBlock(final int index)
        throws InvalidMatrixException {

        // the characteristic polynomial is
        // X^3 - (q0 + q1 + q2) X^2 + (q0 q1 + q0 q2 + q1 q2 - e1^2 - e2^2) X + q0 e2^2 + q2 e1^2 - q0 q1 q2
        final double q0       = main[index];
        final double q1       = main[index + 1];
        final double q2       = main[index + 2];
        final double e12      = squaredSecondary[index];
        final double q1q2Me22 = q1 * q2 - squaredSecondary[index + 1];

        // compute coefficients of the cubic equation as: x^3 + b x^2 + c x + d = 0
        final double b        = -(q0 + q1 + q2);
        final double c        = q0 * q1 + q0 * q2 + q1q2Me22 - e12;
        final double d        = q2 * e12 - q0 * q1q2Me22;

        // solve cubic equation
        final double b2       = b * b;
        final double q        = (3 * c - b2) / 9;
        final double r        = ((9 * c - 2 * b2) * b - 27 * d) / 54;
        final double delta    = q * q * q + r * r;
        if (delta >= 0) {
            // in fact, there are solutions to the equation, but in the context
            // of symmetric realEigenvalues problem, there should be three distinct
            // real roots, so we throw an error if this condition is not met
            throw new InvalidMatrixException("cannot solve degree {0} equation", 3);           
        }
        final double sqrtMq = Math.sqrt(-q);
        final double theta  = Math.acos(r / (-q * sqrtMq));
        final double alpha  = 2 * sqrtMq;
        final double beta   = b / 3;

        double z0 = alpha * Math.cos(theta / 3) - beta;
        double z1 = alpha * Math.cos((theta + 2 * Math.PI) / 3) - beta;
        double z2 = alpha * Math.cos((theta + 4 * Math.PI) / 3) - beta;
        if (z0 < z1) {
            final double t = z0;
            z0 = z1;
            z1 = t;
        }
        if (z1 < z2) {
            final double t = z1;
            z1 = z2;
            z2 = t;
        }
        if (z0 < z1) {
            final double t = z0;
            z0 = z1;
            z1 = t;
        }
        realEigenvalues[index]     = z0;
        realEigenvalues[index + 1] = z1;
        realEigenvalues[index + 2] = z2;

    }

    /**
     * Find realEigenvalues using dqd/dqds algorithms.
     * <p>This implementation is based on Beresford N. Parlett
     * and Osni A. Marques paper <a
     * href="http://www.netlib.org/lapack/lawnspdf/lawn155.pdf">An
     * Implementation of the dqds Algorithm (Positive Case)</a> and on the
     * corresponding LAPACK routine DLASQ2.</p>
     * @param n number of rows of the block
     * @exception InvalidMatrixException if block cannot be diagonalized
     * after 30 * n iterations
     */
    private void processGeneralBlock(final int n)
        throws InvalidMatrixException {

        // check decomposed matrix data range
        double sumOffDiag = 0;
        for (int i = 0; i < n - 1; ++i) {
            final int fourI = 4 * i;
            final double ei = work[fourI + 2];
            sumOffDiag += ei;
        }

        if (sumOffDiag == 0) {
            // matrix is already diagonal
            return;
        }

        // initial checks for splits (see Parlett & Marques section 3.3)
        flipIfWarranted(n, 2);

        // two iterations with Li's test for initial splits
        initialSplits(n);

        // initialize parameters used by goodStep
        tType = 0;
        dMin1 = 0;
        dMin2 = 0;
        dN    = 0;
        dN1   = 0;
        dN2   = 0;
        tau   = 0;

        // process split segments
        int i0 = 0;
        int n0 = n;
        while (n0 > 0) {

            // retrieve shift that was temporarily stored as a negative off-diagonal element
            sigma    = (n0 == n) ? 0 : -work[4 * n0 - 2];
            sigmaLow = 0;

            // find start of a new split segment to process
            double eMin = (i0 == n0) ? 0 : work[4 * n0 - 6];
            double eMax = 0;
            double qMax = work[4 * n0 - 4];
            double qMin = qMax;
            i0 = 0;
            for (int i = 4 * (n0 - 2); i >= 0; i -= 4) {
                if (work[i + 2] <= 0) {
                    i0 = 1 + i / 4;
                    break;
                }
                if (qMin >= 4 * eMax) {
                    qMin = Math.min(qMin, work[i + 4]);
                    eMax = Math.max(eMax, work[i + 2]);
                }
                qMax = Math.max(qMax, work[i] + work[i + 2]);
                eMin = Math.min(eMin, work[i + 2]);
            }
            work[4 * n0 - 2] = eMin;

            // lower bound of Gershgorin disk
            dMin = -Math.max(0, qMin - 2 * Math.sqrt(qMin * eMax));

            pingPong = 0;
            int maxIter = 30 * (n0 - i0);
            for (int k = 0; i0 < n0; ++k) {
                if (k >= maxIter) {
                    throw new InvalidMatrixException(new MaxIterationsExceededException(maxIter));
                }

                // perform one step
                n0 = goodStep(i0, n0);
                pingPong = 1 - pingPong;

                // check for new splits after "ping" steps
                // when the last elements of qd array are very small
                if ((pingPong == 0) && (n0 - i0 > 3) &&
                    (work[4 * n0 - 1] <= TOLERANCE_2 * qMax) &&
                    (work[4 * n0 - 2] <= TOLERANCE_2 * sigma)) {
                    int split = i0 - 1;
                    qMax = work[4 * i0];
                    eMin = work[4 * i0 + 2];
                    double previousEMin = work[4 * i0 + 3];
                    for (int i = 4 * i0; i < 4 * n0 - 11; i += 4) {
                        if ((work[i + 3] <= TOLERANCE_2 * work[i]) &&
                            (work[i + 2] <= TOLERANCE_2 * sigma)) {
                            // insert a split
                            work[i + 2]  = -sigma;
                            split        = i / 4;
                            qMax         = 0;
                            eMin         = work[i + 6];
                            previousEMin = work[i + 7];
                        } else {
                            qMax         = Math.max(qMax, work[i + 4]);
                            eMin         = Math.min(eMin, work[i + 2]);
                            previousEMin = Math.min(previousEMin, work[i + 3]);
                        }
                    }
                    work[4 * n0 - 2] = eMin;
                    work[4 * n0 - 1] = previousEMin;
                    i0 = split + 1;
                }
            }

        }

    }

    /**
     * Perform two iterations with Li's tests for initial splits.
     * @param n number of rows of the matrix to process
     */
    private void initialSplits(final int n) {

        pingPong = 0;
        for (int k = 0; k < 2; ++k) {

            // apply Li's reverse test
            double d = work[4 * (n - 1) + pingPong];
            for (int i = 4 * (n - 2) + pingPong; i >= 0; i -= 4) {
                if (work[i + 2] <= TOLERANCE_2 * d) {
                    work[i + 2] = -0.0;
                    d = work[i];
                } else {
                    d *= work[i] / (d + work[i + 2]);
                }
            }

            // apply dqd plus Li's forward test.
            d = work[pingPong];
            for (int i = 2 + pingPong; i < 4 * n - 2; i += 4) {
                final int j = i - 2 * pingPong - 1;
                work[j] = d + work[i];
                if (work[i] <= TOLERANCE_2 * d) {
                    work[i]     = -0.0;
                    work[j]     = d;
                    work[j + 2] = 0.0;
                    d = work[i + 2];
                } else if ((MathUtils.SAFE_MIN * work[i + 2] < work[j]) &&
                           (MathUtils.SAFE_MIN * work[j] < work[i + 2])) {
                    final double tmp = work[i + 2] / work[j];
                    work[j + 2] = work[i] * tmp;
                    d *= tmp;
                } else {
                    work[j + 2] = work[i + 2] * (work[i] / work[j]);
                    d *= work[i + 2] / work[j];
               }
            }
            work[4 * n - 3 - pingPong] = d;

            // from ping to pong
            pingPong = 1 - pingPong;

        }

    }

    /**
     * Perform one "good" dqd/dqds step.
     * <p>This implementation is based on Beresford N. Parlett
     * and Osni A. Marques paper <a
     * href="http://www.netlib.org/lapack/lawnspdf/lawn155.pdf">An
     * Implementation of the dqds Algorithm (Positive Case)</a> and on the
     * corresponding LAPACK routine DLAZQ3.</p>
     * @param start start index
     * @param end end index
     * @return new end (maybe deflated)
     */
    private int goodStep(final int start, final int end) {

        g = 0.0;

        // step 1: accepting realEigenvalues
        int deflatedEnd = end;
        for (boolean deflating = true; deflating;) {

            if (start >= deflatedEnd) {
                // the array has been completely deflated
                return deflatedEnd;
            }

            final int k = 4 * deflatedEnd + pingPong - 1;

            if ((start == deflatedEnd - 1) ||
                ((start != deflatedEnd - 2) &&
                 ((work[k - 5] <= TOLERANCE_2 * (sigma + work[k - 3])) ||
                  (work[k - 2 * pingPong - 4] <= TOLERANCE_2 * work[k - 7])))) {

                // one eigenvalue found, deflate array
                work[4 * deflatedEnd - 4] = sigma + work[4 * deflatedEnd - 4 + pingPong];
                deflatedEnd -= 1;

            } else if ((start == deflatedEnd - 2) ||
                (work[k - 9] <= TOLERANCE_2 * sigma) ||
                (work[k - 2 * pingPong - 8] <= TOLERANCE_2 * work[k - 11])) {

                // two realEigenvalues found, deflate array
                if (work[k - 3] > work[k - 7]) {
                    final double tmp = work[k - 3];
                    work[k - 3] = work[k - 7];
                    work[k - 7] = tmp;
                }

                if (work[k - 5] > TOLERANCE_2 * work[k - 3]) {
                    double t = 0.5 * ((work[k - 7] - work[k - 3]) + work[k - 5]);
                    double s = work[k - 3] * (work[k - 5] / t);
                    if (s <= t) {
                        s = work[k - 3] * work[k - 5] / (t * (1 + Math.sqrt(1 + s / t)));
                    } else {
                        s = work[k - 3] * work[k - 5] / (t + Math.sqrt(t * (t + s)));                      
                    }
                    t = work[k - 7] + (s + work[k - 5]);
                    work[k - 3] *= work[k - 7] / t;
                    work[k - 7]  = t;
                }
                work[4 * deflatedEnd - 8] = sigma + work[k - 7];
                work[4 * deflatedEnd - 4] = sigma + work[k - 3];
                deflatedEnd -= 2;
            } else {

                // no more realEigenvalues found, we need to iterate
                deflating = false;

            }

        }

        final int l = 4 * deflatedEnd + pingPong - 1;

        // step 2: flip array if needed
        if ((dMin <= 0) || (deflatedEnd < end)) {
            if (flipIfWarranted(deflatedEnd, 1)) {
                dMin2 = Math.min(dMin2, work[l - 1]);
                work[l - 1] =
                    Math.min(work[l - 1],
                             Math.min(work[3 + pingPong], work[7 + pingPong]));
                work[l - 2 * pingPong] =
                    Math.min(work[l - 2 * pingPong],
                             Math.min(work[6 + pingPong], work[6 + pingPong]));
                qMax  = Math.max(qMax, Math.max(work[3 + pingPong], work[7 + pingPong]));
                dMin  = -0.0;
            }
        }

        if ((dMin < 0) ||
            (MathUtils.SAFE_MIN * qMax < Math.min(work[l - 1],
                                                  Math.min(work[l - 9],
                                                           dMin2 + work[l - 2 * pingPong])))) {
            // step 3: choose a shift
            computeShiftIncrement(start, deflatedEnd, end - deflatedEnd);

            // step 4a: dqds
            for (boolean loop = true; loop;) {

                // perform one dqds step with the chosen shift
                dqds(start, deflatedEnd);

                // check result of the dqds step
                if ((dMin >= 0) && (dMin1 > 0)) {
                    // the shift was good
                    updateSigma(tau);
                    return deflatedEnd;
                } else if ((dMin < 0.0) &&
                           (dMin1 > 0.0) &&
                           (work[4 * deflatedEnd - 5 - pingPong] < TOLERANCE * (sigma + dN1)) &&
                           (Math.abs(dN) < TOLERANCE * sigma)) {
                   // convergence hidden by negative DN.
                    work[4 * deflatedEnd - 3 - pingPong] = 0.0;
                    dMin = 0.0;
                    updateSigma(tau);
                    return deflatedEnd;
                } else if (dMin < 0.0) {
                    // tau too big. Select new tau and try again.
                    if (tType < -22) {
                        // failed twice. Play it safe.
                        tau = 0.0;
                    } else if (dMin1 > 0.0) {
                        // late failure. Gives excellent shift.
                        tau = (tau + dMin) * (1.0 - 2.0 * MathUtils.EPSILON);
                        tType -= 11;
                    } else {
                        // early failure. Divide by 4.
                        tau *= 0.25;
                        tType -= 12;
                    }
                } else if (Double.isNaN(dMin)) {
                    tau = 0.0;
                } else {
                    // possible underflow. Play it safe.
                    loop = false;
                }
            }

        }

        // perform a dqd step (i.e. no shift)
        dqd(start, deflatedEnd);

        return deflatedEnd;

    }

    /**
     * Flip qd array if warranted.
     * @param n number of rows in the block
     * @param step within the array (1 for flipping all elements, 2 for flipping
     * only every other element)
     * @return true if qd array was flipped
     */
    private boolean flipIfWarranted(final int n, final int step) {
        if (1.5 * work[pingPong] < work[4 * (n - 1) + pingPong]) {
            // flip array
            for (int i = 0, j = 4 * n - 1; i < j; i += 4, j -= 4) {
                for (int k = 0; k < 4; k += step) {
                    final double tmp = work[i + k];
                    work[i + k] = work[j - k];
                    work[j - k] = tmp;
                }
            }
            return true;
        }
        return false;
    }

    /**
     * Compute an interval containing all realEigenvalues of a block.
     * @param index index of the first row of the block
     * @param n number of rows of the block
     * @return an interval containing the realEigenvalues
     */
    private double[] eigenvaluesRange(final int index, final int n) {

        // find the bounds of the spectra of the local block
        final int lowerStart = 4 * main.length;
        final int upperStart = 5 * main.length;
        double lower = Double.POSITIVE_INFINITY;
        double upper = Double.NEGATIVE_INFINITY;
        for (int i = 0; i < n; ++i) {
            lower = Math.min(lower, work[lowerStart + index +i]);
            upper = Math.max(upper, work[upperStart + index +i]);
        }

        // set thresholds
        final double tNorm = Math.max(Math.abs(lower), Math.abs(upper));
        final double relativeTolerance = Math.sqrt(MathUtils.EPSILON);
        final double absoluteTolerance = 4 * minPivot;
        final int maxIter =
            2 + (int) ((Math.log(tNorm + minPivot) - Math.log(minPivot)) / Math.log(2.0));
        final double margin = 2 * (tNorm * MathUtils.EPSILON * n + 2 * minPivot);

        // search lower eigenvalue
        double left  = lower - margin;
        double right = upper + margin;
        for (int i = 0; i < maxIter; ++i) {

            final double range = right - left;
            if ((range < absoluteTolerance) ||
                (range < relativeTolerance * Math.max(Math.abs(left), Math.abs(right)))) {
                // search has converged
                break;
            }

            final double middle = 0.5 * (left + right);
            if (countEigenValues(middle, index, n) >= 1) {
                right = middle;
            } else {
                left = middle;
            }

        }
        lower = Math.max(lower, left - 100 * MathUtils.EPSILON * Math.abs(left));

        // search upper eigenvalue
        left  = lower - margin;
        right = upper + margin;
        for (int i = 0; i < maxIter; ++i) {

            final double range = right - left;
            if ((range < absoluteTolerance) ||
                (range < relativeTolerance * Math.max(Math.abs(left), Math.abs(right)))) {
                // search has converged
                break;
            }

            final double middle = 0.5 * (left + right);
            if (countEigenValues(middle, index, n) >= n) {
                right = middle;
            } else {
                left = middle;
            }

        }
        upper = Math.min(upper, right + 100 * MathUtils.EPSILON * Math.abs(right));

        return new double[] { lower, upper };

    }

    /**
     * Count the number of realEigenvalues below a point.
     * @param t value below which we must count the number of realEigenvalues
     * @param index index of the first row of the block
     * @param n number of rows of the block
     * @return number of realEigenvalues smaller than t
     */
    private int countEigenValues(final double t, final int index, final int n) {
        double ratio = main[index] - t;
        int count = (ratio > 0) ? 0 : 1;
        for (int i = 1; i < n; ++i) {
            ratio = main[index + i] - squaredSecondary[index + i - 1] / ratio - t;
            if (ratio <= 0) {
                ++count;
            }
        }
        return count;
    }

    /**
     * Decompose the shifted tridiagonal matrix T-&lambda;I as LDL<sup>T</sup>.
     * <p>A shifted symmetric tridiagonal matrix T can be decomposed as
     * LDL<sup>T</sup> where L is a lower bidiagonal matrix with unit diagonal
     * and D is a diagonal matrix. This method is an implementation of
     * algorithm 4.4.7 from Dhillon's thesis.</p>
     * @param lambda shift to add to the matrix before decomposing it
     * to ensure it is positive definite
     * @param index index of the first row of the block
     * @param n number of rows of the block
     */
    private void ldlTDecomposition(final double lambda, final int index, final int n) {
        double di = main[index] - lambda;
        work[0] = Math.abs(di);
        for (int i = 1; i < n; ++i) {
            final int    fourI = 4 * i;
            final double eiM1  = secondary[index + i - 1];
            final double ratio = eiM1 / di;
            work[fourI - 2] = ratio * ratio * Math.abs(di);
            di = (main[index + i] - lambda) - eiM1 * ratio;
            work[fourI] = Math.abs(di);
        }
    }

    /**
     * Perform a dqds step, using current shift increment.
     * <p>This implementation is a translation of the LAPACK routine DLASQ5.</p>
     * @param start start index
     * @param end end index
     */
    private void dqds(final int start, final int end) {

        eMin = work[4 * start + pingPong + 4];
        double d = work[4 * start + pingPong] - tau;
        dMin = d;
        dMin1 = -work[4 * start + pingPong];

        if (pingPong == 0) {
            for (int j4 = 4 * start + 3; j4 <= 4 * (end - 3); j4 += 4) {
                work[j4 - 2] = d + work[j4 - 1];
                final double tmp = work[j4 + 1] / work[j4 - 2];
                d = d * tmp - tau;
                dMin = Math.min(dMin, d);
                work[j4] = work[j4 - 1] * tmp;
                eMin = Math.min(work[j4], eMin);
            }
        } else {
            for (int j4 = 4 * start + 3; j4 <= 4 * (end - 3); j4 += 4) {
                work[j4 - 3] = d + work[j4];
                final double tmp = work[j4 + 2] / work[j4 - 3];
                d = d * tmp - tau;
                dMin = Math.min(dMin, d);
                work[j4 - 1] = work[j4] * tmp;
                eMin = Math.min(work[j4 - 1], eMin);
            }
        }

        // unroll last two steps.
        dN2 = d;
        dMin2 = dMin;
        int j4 = 4 * (end - 2) - pingPong - 1;
        int j4p2 = j4 + 2 * pingPong - 1;
        work[j4 - 2] = dN2 + work[j4p2];
        work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
        dN1 = work[j4p2 + 2] * (dN2 / work[j4 - 2]) - tau;
        dMin = Math.min(dMin, dN1);

        dMin1 = dMin;
        j4 = j4 + 4;
        j4p2 = j4 + 2 * pingPong - 1;
        work[j4 - 2] = dN1 + work[j4p2];
        work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
        dN = work[j4p2 + 2] * (dN1 / work[j4 - 2]) - tau;
        dMin = Math.min(dMin, dN);

        work[j4 + 2] = dN;
        work[4 * end - pingPong - 1] = eMin;

    }


    /**
     * Perform a dqd step.
     * <p>This implementation is a translation of the LAPACK routine DLASQ6.</p>
     * @param start start index
     * @param end end index
     */
    private void dqd(final int start, final int end) {

        eMin = work[4 * start + pingPong + 4];
        double d = work[4 * start + pingPong];
        dMin = d;

        if (pingPong == 0) {
            for (int j4 = 4 * start + 3; j4 < 4 * (end - 3); j4 += 4) {
                work[j4 - 2] = d + work[j4 - 1];
                if (work[j4 - 2] == 0.0) {
                    work[j4] = 0.0;
                    d = work[j4 + 1];
                    dMin = d;
                    eMin = 0.0;
                } else if ((MathUtils.SAFE_MIN * work[j4 + 1] < work[j4 - 2]) &&
                           (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4 + 1])) {
                    final double tmp = work[j4 + 1] / work[j4 - 2];
                    work[j4] = work[j4 - 1] * tmp;
                    d *= tmp;
                } else {
                    work[j4] = work[j4 + 1] * (work[j4 - 1] / work[j4 - 2]);
                    d *= work[j4 + 1] / work[j4 - 2];
                }
                dMin = Math.min(dMin, d);
                eMin = Math.min(eMin, work[j4]);
            }
        } else {
            for (int j4 = 4 * start + 3; j4 < 4 * (end - 3); j4 += 4) {
                work[j4 - 3] = d + work[j4];
                if (work[j4 - 3] == 0.0) {
                    work[j4 - 1] = 0.0;
                    d = work[j4 + 2];
                    dMin = d;
                    eMin = 0.0;
                } else if ((MathUtils.SAFE_MIN * work[j4 + 2] < work[j4 - 3]) &&
                           (MathUtils.SAFE_MIN * work[j4 - 3] < work[j4 + 2])) {
                    final double tmp = work[j4 + 2] / work[j4 - 3];
                    work[j4 - 1] = work[j4] * tmp;
                    d *= tmp;
                } else {
                    work[j4 - 1] = work[j4 + 2] * (work[j4] / work[j4 - 3]);
                    d *= work[j4 + 2] / work[j4 - 3];
                }
                dMin = Math.min(dMin, d);
                eMin = Math.min(eMin, work[j4 - 1]);
            }
        }

        // Unroll last two steps
        dN2   = d;
        dMin2 = dMin;
        int j4 = 4 * (end - 2) - pingPong - 1;
        int j4p2 = j4 + 2 * pingPong - 1;
        work[j4 - 2] = dN2 + work[j4p2];
        if (work[j4 - 2] == 0.0) {
            work[j4] = 0.0;
            dN1  = work[j4p2 + 2];
            dMin = dN1;
            eMin = 0.0;
        } else if ((MathUtils.SAFE_MIN * work[j4p2 + 2] < work[j4 - 2]) &&
                   (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4p2 + 2])) {
            final double tmp = work[j4p2 + 2] / work[j4 - 2];
            work[j4] = work[j4p2] * tmp;
            dN1 = dN2 * tmp;
        } else {
            work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
            dN1 = work[j4p2 + 2] * (dN2 / work[j4 - 2]);
        }
        dMin = Math.min(dMin, dN1);

        dMin1 = dMin;
        j4 = j4 + 4;
        j4p2 = j4 + 2 * pingPong - 1;
        work[j4 - 2] = dN1 + work[j4p2];
        if (work[j4 - 2] == 0.0) {
            work[j4] = 0.0;
            dN   = work[j4p2 + 2];
            dMin = dN;
            eMin = 0.0;
        } else if ((MathUtils.SAFE_MIN * work[j4p2 + 2] < work[j4 - 2]) &&
                   (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4p2 + 2])) {
            final double tmp = work[j4p2 + 2] / work[j4 - 2];
            work[j4] = work[j4p2] * tmp;
            dN = dN1 * tmp;
        } else {
            work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
            dN = work[j4p2 + 2] * (dN1 / work[j4 - 2]);
        }
        dMin = Math.min(dMin, dN);

        work[j4 + 2] = dN;
        work[4 * end - pingPong - 1] = eMin;

    }

    /**
     * Compute the shift increment as an estimate of the smallest eigenvalue.
     * <p>This implementation is a translation of the LAPACK routine DLAZQ4.</p>
     * @param start start index
     * @param end end index
     * @param deflated number of realEigenvalues just deflated
     */
    private void computeShiftIncrement(final int start, final int end, final int deflated) {

        final double cnst1 = 0.563;
        final double cnst2 = 1.010;
        final double cnst3 = 1.05;

        // a negative dMin forces the shift to take that absolute value
        // tType records the type of shift.
        if (dMin <= 0.0) {
            tau = -dMin;
            tType = -1;
            return;
        }

        int nn = 4 * end + pingPong - 1;
        switch (deflated) {

        case 0 : // no realEigenvalues deflated. 
            if (dMin == dN || dMin == dN1) {

                double b1 = Math.sqrt(work[nn - 3]) * Math.sqrt(work[nn - 5]);
                double b2 = Math.sqrt(work[nn - 7]) * Math.sqrt(work[nn - 9]);
                double a2 = work[nn - 7] + work[nn - 5];

                if (dMin == dN && dMin1 == dN1) {
                    // cases 2 and 3. 
                    final double gap2 = dMin2 - a2 - dMin2 * 0.25;
                    final double gap1 = a2 - dN - ((gap2 > 0.0 && gap2 > b2) ? (b2 / gap2) * b2 : (b1 + b2));
                    if (gap1 > 0.0 && gap1 > b1) {
                        tau   = Math.max(dN - (b1 / gap1) * b1, 0.5 * dMin);
                        tType = -2;
                    } else {
                        double s = 0.0;
                        if (dN > b1) {
                            s = dN - b1;
                        }
                        if (a2 > (b1 + b2)) {
                            s = Math.min(s, a2 - (b1 + b2));
                        }
                        tau   = Math.max(s, 0.333 * dMin);
                        tType = -3;
                    }
                } else {
                    // case 4.
                    tType = -4;
                    double s = 0.25 * dMin;
                    double gam;
                    int np;
                    if (dMin == dN) {
                        gam = dN;
                        a2 = 0.0;
                        if (work[nn - 5]  >  work[nn - 7]) {
                            return;
                        }
                        b2 = work[nn - 5] / work[nn - 7];
                        np = nn - 9;
                    } else {
                        np = nn - 2 * pingPong;
                        b2 = work[np - 2];
                        gam = dN1;
                        if (work[np - 4]  >  work[np - 2]) {
                            return;
                        }
                        a2 = work[np - 4] / work[np - 2];
                        if (work[nn - 9]  >  work[nn - 11]) {
                            return;
                        }
                        b2 = work[nn - 9] / work[nn - 11];
                        np = nn - 13;
                    }

                    // approximate contribution to norm squared from i < nn-1.
                    a2 = a2 + b2;
                    for (int i4 = np; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
                        if(b2 == 0.0) {
                            break;
                        }
                        b1 = b2;
                        if (work[i4]  >  work[i4 - 2]) {
                            return;
                        }
                        b2 = b2 * (work[i4] / work[i4 - 2]);
                        a2 = a2 + b2;
                        if (100 * Math.max(b2, b1) < a2 || cnst1 < a2) {
                            break;
                        }
                    }
                    a2 = cnst3 * a2;

                    // rayleigh quotient residual bound.
                    if (a2 < cnst1) {
                        s = gam * (1 - Math.sqrt(a2)) / (1 + a2);
                    }
                    tau = s;

                }
            } else if (dMin == dN2) {

                // case 5.
                tType = -5;
                double s = 0.25 * dMin;

                // compute contribution to norm squared from i > nn-2.
                final int np = nn - 2 * pingPong;
                double b1 = work[np - 2];
                double b2 = work[np - 6];
                final double gam = dN2;
                if (work[np - 8] > b2 || work[np - 4] > b1) {
                    return;
                }
                double a2 = (work[np - 8] / b2) * (1 + work[np - 4] / b1);

                // approximate contribution to norm squared from i < nn-2.
                if (end - start > 2) {
                    b2 = work[nn - 13] / work[nn - 15];
                    a2 = a2 + b2;
                    for (int i4 = nn - 17; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
                        if (b2 == 0.0) {
                            break;
                        }
                        b1 = b2;
                        if (work[i4]  >  work[i4 - 2]) {
                            return;
                        }
                        b2 = b2 * (work[i4] / work[i4 - 2]);
                        a2 = a2 + b2;
                        if (100 * Math.max(b2, b1) < a2 || cnst1 < a2)  {
                            break;
                        }
                    }
                    a2 = cnst3 * a2;
                }

                if (a2 < cnst1) {
                    tau = gam * (1 - Math.sqrt(a2)) / (1 + a2);
                } else {
                    tau = s;
                }

            } else {

                // case 6, no information to guide us.
                if (tType == -6) {
                    g += 0.333 * (1 - g);
                } else if (tType == -18) {
                    g = 0.25 * 0.333;
                } else {
                    g = 0.25;
                }
                tau   = g * dMin;
                tType = -6;

            }
            break;

        case 1 : // one eigenvalue just deflated. use dMin1, dN1 for dMin and dN.
            if (dMin1 == dN1 && dMin2 == dN2) { 

                // cases 7 and 8.
                tType = -7;
                double s = 0.333 * dMin1;
                if (work[nn - 5] > work[nn - 7]) {
                    return;
                }
                double b1 = work[nn - 5] / work[nn - 7];
                double b2 = b1;
                if (b2 != 0.0) {
                    for (int i4 = 4 * end - 10 + pingPong; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
                        final double oldB1 = b1;
                        if (work[i4] > work[i4 - 2]) {
                            return;
                        }
                        b1 = b1 * (work[i4] / work[i4 - 2]);
                        b2 = b2 + b1;
                        if (100 * Math.max(b1, oldB1) < b2) {
                            break;
                        }
                    }
                }
                b2 = Math.sqrt(cnst3 * b2);
                final double a2 = dMin1 / (1 + b2 * b2);
                final double gap2 = 0.5 * dMin2 - a2;
                if (gap2 > 0.0 && gap2 > b2 * a2) {
                    tau = Math.max(s, a2 * (1 - cnst2 * a2 * (b2 / gap2) * b2));
                } else {
                    tau = Math.max(s, a2 * (1 - cnst2 * b2));
                    tType = -8;
                }
            } else {

                // case 9.
                tau = 0.25 * dMin1;
                if (dMin1 == dN1) {
                    tau = 0.5 * dMin1;
                }
                tType = -9;
            }
            break;

        case 2 : // two realEigenvalues deflated. use dMin2, dN2 for dMin and dN.

            // cases 10 and 11.
            if (dMin2 == dN2 && 2 * work[nn - 5] < work[nn - 7]) { 
                tType = -10;
                final double s = 0.333 * dMin2;
                if (work[nn - 5] > work[nn - 7]) {
                    return;
                }
                double b1 = work[nn - 5] / work[nn - 7];
                double b2 = b1;
                if (b2 != 0.0){
                    for (int i4 = 4 * end - 9 + pingPong; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
                        if (work[i4] > work[i4 - 2]) {
                            return;
                        }
                        b1 *= work[i4] / work[i4 - 2];
                        b2 += b1;
                        if (100 * b1 < b2) {
                            break;
                        }
                    }
                }
                b2 = Math.sqrt(cnst3 * b2);
                final double a2 = dMin2 / (1 + b2 * b2);
                final double gap2 = work[nn - 7] + work[nn - 9] -
                Math.sqrt(work[nn - 11]) * Math.sqrt(work[nn - 9]) - a2;
                if (gap2 > 0.0 && gap2 > b2 * a2) {
                    tau = Math.max(s, a2 * (1 - cnst2 * a2 * (b2 / gap2) * b2));
                } else {
                    tau = Math.max(s, a2 * (1 - cnst2 * b2));
                }
            } else {
                tau   = 0.25 * dMin2;
                tType = -11;
            }
            break;

        default : // case 12, more than two realEigenvalues deflated. no information.
            tau   = 0.0;
            tType = -12;
        }

    }

    /**
     * Update sigma.
     * @param tau shift to apply to sigma
     */
    private void updateSigma(final double tau) {
        // BEWARE: do NOT attempt to simplify the following statements
        // the expressions below take care to accumulate the part of sigma
        // that does not fit within a double variable into sigmaLow
        if (tau < sigma) {
            sigmaLow += tau;
            final double t = sigma + sigmaLow;
            sigmaLow -= t - sigma;
            sigma = t;
        } else {
            final double t = sigma + tau;
            sigmaLow += sigma - (t - tau);
            sigma = t;
        }
    }

    /**
     * Find eigenvectors.
     */
    private void findEigenVectors() {

        final int m = main.length;
        eigenvectors = new ArrayRealVector[m];

        // perform an initial non-shifted LDLt decomposition
        final double[] d = new double[m];
        final double[] l = new double[m - 1];
        double di = main[0];
        d[0] = di;
        for (int i = 1; i < m; ++i) {
            final double eiM1  = secondary[i - 1];
            final double ratio = eiM1 / di;
            di       = main[i] - eiM1 * ratio;
            l[i - 1] = ratio;
            d[i]     = di;
        }

        // compute eigenvectors
        for (int i = 0; i < m; ++i) {
            eigenvectors[i] = findEigenvector(realEigenvalues[i], d, l);
        }

    }

    /**
     * Find an eigenvector corresponding to an eigenvalue, using bidiagonals.
     * <p>This method corresponds to algorithm X from Dhillon's thesis.</p>
     * 
     * @param eigenvalue eigenvalue for which eigenvector is desired
     * @param d diagonal elements of the initial non-shifted D matrix
     * @param l off-diagonal elements of the initial non-shifted L matrix
     * @return an eigenvector
     */
    private ArrayRealVector findEigenvector(final double eigenvalue,
                                           final double[] d, final double[] l) {

        // compute the LDLt and UDUt decompositions of the
        // perfectly shifted tridiagonal matrix
        final int m = main.length;
        stationaryQuotientDifferenceWithShift(d, l, eigenvalue);
        progressiveQuotientDifferenceWithShift(d, l, eigenvalue);

        // select the twist index leading to
        // the least diagonal element in the twisted factorization
        int r = m - 1;
        double minG = Math.abs(work[6 * r] + work[6 * r + 3] + eigenvalue);
        for (int i = 0, sixI = 0; i < m - 1; ++i, sixI += 6) {
            final double g = work[sixI] + d[i] * work[sixI + 9] / work[sixI + 10];
            final double absG = Math.abs(g);
            if (absG < minG) {
                r = i;
                minG = absG;
            }
        }

        // solve the singular system by ignoring the equation
        // at twist index and propagating upwards and downwards
        double[] eigenvector = new double[m];
        double n2 = 1;
        eigenvector[r] = 1;
        double z = 1;
        for (int i = r - 1; i >= 0; --i) {
            z *= -work[6 * i + 2];
            eigenvector[i] = z;
            n2 += z * z;
        }
        z = 1;
        for (int i = r + 1; i < m; ++i) {
            z *= -work[6 * i - 1];
            eigenvector[i] = z;
            n2 += z * z;
        }

        // normalize vector
        final double inv = 1.0 / Math.sqrt(n2);
        for (int i = 0; i < m; ++i) {
            eigenvector[i] *= inv;
        }

        return (transformer == null) ?
               new ArrayRealVector(eigenvector, false) :
               new ArrayRealVector(transformer.getQ().operate(eigenvector), false);

    }

    /**
     * Decompose matrix LDL<sup>T</sup> - &lambda; I as
     * L<sub>+</sub>D<sub>+</sub>L<sub>+</sub><sup>T</sup>.
     * <p>This method corresponds to algorithm 4.4.3 (dstqds) from Dhillon's thesis.</p>
     * @param d diagonal elements of D,
     * @param l off-diagonal elements of L
     * @param lambda shift to apply
     */
    private void stationaryQuotientDifferenceWithShift(final double[] d, final double[] l,
                                                       final double lambda) {
        final int nM1 = d.length - 1;
        double si = -lambda;
        for (int i = 0, sixI = 0; i < nM1; ++i, sixI += 6) {
            final double di   = d[i];
            final double li   = l[i];
            final double diP1 = di + si;
            final double liP1 = li * di / diP1;
            work[sixI]        = si;
            work[sixI + 1]    = diP1;
            work[sixI + 2]    = liP1;
            si = li * liP1 * si - lambda;
        }
        work[6 * nM1 + 1] = d[nM1] + si;
        work[6 * nM1]     = si;
    }

    /**
     * Decompose matrix LDL<sup>T</sup> - &lambda; I as
     * U<sub>-</sub>D<sub>-</sub>U<sub>-</sub><sup>T</sup>.
     * <p>This method corresponds to algorithm 4.4.5 (dqds) from Dhillon's thesis.</p>
     * @param d diagonal elements of D
     * @param l off-diagonal elements of L
     * @param lambda shift to apply
     */
    private void progressiveQuotientDifferenceWithShift(final double[] d, final double[] l,
                                                        final double lambda) {
        final int nM1 = d.length - 1;
        double pi = d[nM1] - lambda;
        for (int i = nM1 - 1, sixI = 6 * i; i >= 0; --i, sixI -= 6) {
            final double di   = d[i];
            final double li   = l[i];
            final double diP1 = di * li * li + pi;
            final double t    = di / diP1;
            work[sixI +  9]   = pi;
            work[sixI + 10]   = diP1;
            work[sixI +  5]   = li * t;
            pi = pi * t - lambda;
        }
        work[3] = pi;
        work[4] = pi;
    }

}
