当前路径：Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php
<?php
/**
 *    @package JAMA
 *
 *    Class to obtain eigenvalues and eigenvectors of a real matrix.
 *
 *    If A is symmetric, then A = V*D*V' where the eigenvalue matrix D
 *    is diagonal and the eigenvector matrix V is orthogonal (i.e.
 *    A = V.times(D.times(V.transpose())) and V.times(V.transpose())
 *    equals the identity matrix).
 *
 *    If A is not symmetric, then the eigenvalue matrix D is block diagonal
 *    with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
 *    lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda].  The
 *    columns of V represent the eigenvectors in the sense that A*V = V*D,
 *    i.e. A.times(V) equals V.times(D).  The matrix V may be badly
 *    conditioned, or even singular, so the validity of the equation
 *    A = V*D*inverse(V) depends upon V.cond().
 *
 *    @author  Paul Meagher
 *    @license PHP v3.0
 *    @version 1.1
 */
class EigenvalueDecomposition
{
    /**
     *    Row and column dimension (square matrix).
     *    @var int
     */
    private $n;

    /**
     *    Internal symmetry flag.
     *    @var int
     */
    private $issymmetric;

    /**
     *    Arrays for internal storage of eigenvalues.
     *    @var array
     */
    private $d = array();
    private $e = array();

    /**
     *    Array for internal storage of eigenvectors.
     *    @var array
     */
    private $V = array();

    /**
    *    Array for internal storage of nonsymmetric Hessenberg form.
    *    @var array
    */
    private $H = array();

    /**
    *    Working storage for nonsymmetric algorithm.
    *    @var array
    */
    private $ort;

    /**
    *    Used for complex scalar division.
    *    @var float
    */
    private $cdivr;
    private $cdivi;

    /**
     *    Symmetric Householder reduction to tridiagonal form.
     *
     *    @access private
     */
    private function tred2()
    {
        //  This is derived from the Algol procedures tred2 by
        //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
        //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
        //  Fortran subroutine in EISPACK.
        $this->d = $this->V[$this->n-1];
        // Householder reduction to tridiagonal form.
        for ($i = $this->n-1; $i > 0; --$i) {
            $i_ = $i -1;
            // Scale to avoid under/overflow.
            $h = $scale = 0.0;
            $scale += array_sum(array_map(abs, $this->d));
            if ($scale == 0.0) {
                $this->e[$i] = $this->d[$i_];
                $this->d = array_slice($this->V[$i_], 0, $i_);
                for ($j = 0; $j < $i; ++$j) {
                    $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
                }
            } else {
                // Generate Householder vector.
                for ($k = 0; $k < $i; ++$k) {
                    $this->d[$k] /= $scale;
                    $h += pow($this->d[$k], 2);
                }
                $f = $this->d[$i_];
                $g = sqrt($h);
                if ($f > 0) {
                    $g = -$g;
                }
                $this->e[$i] = $scale * $g;
                $h = $h - $f * $g;
                $this->d[$i_] = $f - $g;
                for ($j = 0; $j < $i; ++$j) {
                    $this->e[$j] = 0.0;
                }
                // Apply similarity transformation to remaining columns.
                for ($j = 0; $j < $i; ++$j) {
                    $f = $this->d[$j];
                    $this->V[$j][$i] = $f;
                    $g = $this->e[$j] + $this->V[$j][$j] * $f;
                    for ($k = $j+1; $k <= $i_; ++$k) {
                        $g += $this->V[$k][$j] * $this->d[$k];
                        $this->e[$k] += $this->V[$k][$j] * $f;
                    }
                    $this->e[$j] = $g;
                }
                $f = 0.0;
                for ($j = 0; $j < $i; ++$j) {
                    $this->e[$j] /= $h;
                    $f += $this->e[$j] * $this->d[$j];
                }
                $hh = $f / (2 * $h);
                for ($j=0; $j < $i; ++$j) {
                    $this->e[$j] -= $hh * $this->d[$j];
                }
                for ($j = 0; $j < $i; ++$j) {
                    $f = $this->d[$j];
                    $g = $this->e[$j];
                    for ($k = $j; $k <= $i_; ++$k) {
                        $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
                    }
                    $this->d[$j] = $this->V[$i-1][$j];
                    $this->V[$i][$j] = 0.0;
                }
            }
            $this->d[$i] = $h;
        }

        // Accumulate transformations.
        for ($i = 0; $i < $this->n-1; ++$i) {
            $this->V[$this->n-1][$i] = $this->V[$i][$i];
            $this->V[$i][$i] = 1.0;
            $h = $this->d[$i+1];
            if ($h != 0.0) {
                for ($k = 0; $k <= $i; ++$k) {
                    $this->d[$k] = $this->V[$k][$i+1] / $h;
                }
                for ($j = 0; $j <= $i; ++$j) {
                    $g = 0.0;
                    for ($k = 0; $k <= $i; ++$k) {
                        $g += $this->V[$k][$i+1] * $this->V[$k][$j];
                    }
                    for ($k = 0; $k <= $i; ++$k) {
                        $this->V[$k][$j] -= $g * $this->d[$k];
                    }
                }
            }
            for ($k = 0; $k <= $i; ++$k) {
                $this->V[$k][$i+1] = 0.0;
            }
        }

        $this->d = $this->V[$this->n-1];
        $this->V[$this->n-1] = array_fill(0, $j, 0.0);
        $this->V[$this->n-1][$this->n-1] = 1.0;
        $this->e[0] = 0.0;
    }

    /**
     *    Symmetric tridiagonal QL algorithm.
     *
     *    This is derived from the Algol procedures tql2, by
     *    Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
     *    Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
     *    Fortran subroutine in EISPACK.
     *
     *    @access private
     */
    private function tql2()
    {
        for ($i = 1; $i < $this->n; ++$i) {
            $this->e[$i-1] = $this->e[$i];
        }
        $this->e[$this->n-1] = 0.0;
        $f = 0.0;
        $tst1 = 0.0;
        $eps  = pow(2.0, -52.0);

        for ($l = 0; $l < $this->n; ++$l) {
            // Find small subdiagonal element
            $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
            $m = $l;
            while ($m < $this->n) {
                if ...