EigenvalueDecomposition.php in Loft Data Grids 7.2

File

vendor/phpoffice/phpexcel/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 (abs($this->e[$m]) <= $eps * $tst1) {
          break;
        }
        ++$m;
      }

      // If m == l, $this->d[l] is an eigenvalue,
      // otherwise, iterate.
      if ($m > $l) {
        $iter = 0;
        do {

          // Could check iteration count here.
          $iter += 1;

          // Compute implicit shift
          $g = $this->d[$l];
          $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]);
          $r = hypo($p, 1.0);
          if ($p < 0) {
            $r *= -1;
          }
          $this->d[$l] = $this->e[$l] / ($p + $r);
          $this->d[$l + 1] = $this->e[$l] * ($p + $r);
          $dl1 = $this->d[$l + 1];
          $h = $g - $this->d[$l];
          for ($i = $l + 2; $i < $this->n; ++$i) {
            $this->d[$i] -= $h;
          }
          $f += $h;

          // Implicit QL transformation.
          $p = $this->d[$m];
          $c = 1.0;
          $c2 = $c3 = $c;
          $el1 = $this->e[$l + 1];
          $s = $s2 = 0.0;
          for ($i = $m - 1; $i >= $l; --$i) {
            $c3 = $c2;
            $c2 = $c;
            $s2 = $s;
            $g = $c * $this->e[$i];
            $h = $c * $p;
            $r = hypo($p, $this->e[$i]);
            $this->e[$i + 1] = $s * $r;
            $s = $this->e[$i] / $r;
            $c = $p / $r;
            $p = $c * $this->d[$i] - $s * $g;
            $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]);

            // Accumulate transformation.
            for ($k = 0; $k < $this->n; ++$k) {
              $h = $this->V[$k][$i + 1];
              $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h;
              $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
            }
          }
          $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
          $this->e[$l] = $s * $p;
          $this->d[$l] = $c * $p;

          // Check for convergence.
        } while (abs($this->e[$l]) > $eps * $tst1);
      }
      $this->d[$l] = $this->d[$l] + $f;
      $this->e[$l] = 0.0;
    }

    // Sort eigenvalues and corresponding vectors.
    for ($i = 0; $i < $this->n - 1; ++$i) {
      $k = $i;
      $p = $this->d[$i];
      for ($j = $i + 1; $j < $this->n; ++$j) {
        if ($this->d[$j] < $p) {
          $k = $j;
          $p = $this->d[$j];
        }
      }
      if ($k != $i) {
        $this->d[$k] = $this->d[$i];
        $this->d[$i] = $p;
        for ($j = 0; $j < $this->n; ++$j) {
          $p = $this->V[$j][$i];
          $this->V[$j][$i] = $this->V[$j][$k];
          $this->V[$j][$k] = $p;
        }
      }
    }
  }

  /**
   *	Nonsymmetric reduction to Hessenberg form.
   *
   *	This is derived from the Algol procedures orthes and ortran,
   *	by Martin and Wilkinson, Handbook for Auto. Comp.,
   *	Vol.ii-Linear Algebra, and the corresponding
   *	Fortran subroutines in EISPACK.
   *
   *	@access private
   */
  private function orthes() {
    $low = 0;
    $high = $this->n - 1;
    for ($m = $low + 1; $m <= $high - 1; ++$m) {

      // Scale column.
      $scale = 0.0;
      for ($i = $m; $i <= $high; ++$i) {
        $scale = $scale + abs($this->H[$i][$m - 1]);
      }
      if ($scale != 0.0) {

        // Compute Householder transformation.
        $h = 0.0;
        for ($i = $high; $i >= $m; --$i) {
          $this->ort[$i] = $this->H[$i][$m - 1] / $scale;
          $h += $this->ort[$i] * $this->ort[$i];
        }
        $g = sqrt($h);
        if ($this->ort[$m] > 0) {
          $g *= -1;
        }
        $h -= $this->ort[$m] * $g;
        $this->ort[$m] -= $g;

        // Apply Householder similarity transformation
        // H = (I -u * u' / h) * H * (I -u * u') / h)
        for ($j = $m; $j < $this->n; ++$j) {
          $f = 0.0;
          for ($i = $high; $i >= $m; --$i) {
            $f += $this->ort[$i] * $this->H[$i][$j];
          }
          $f /= $h;
          for ($i = $m; $i <= $high; ++$i) {
            $this->H[$i][$j] -= $f * $this->ort[$i];
          }
        }
        for ($i = 0; $i <= $high; ++$i) {
          $f = 0.0;
          for ($j = $high; $j >= $m; --$j) {
            $f += $this->ort[$j] * $this->H[$i][$j];
          }
          $f = $f / $h;
          for ($j = $m; $j <= $high; ++$j) {
            $this->H[$i][$j] -= $f * $this->ort[$j];
          }
        }
        $this->ort[$m] = $scale * $this->ort[$m];
        $this->H[$m][$m - 1] = $scale * $g;
      }
    }

    // Accumulate transformations (Algol's ortran).
    for ($i = 0; $i < $this->n; ++$i) {
      for ($j = 0; $j < $this->n; ++$j) {
        $this->V[$i][$j] = $i == $j ? 1.0 : 0.0;
      }
    }
    for ($m = $high - 1; $m >= $low + 1; --$m) {
      if ($this->H[$m][$m - 1] != 0.0) {
        for ($i = $m + 1; $i <= $high; ++$i) {
          $this->ort[$i] = $this->H[$i][$m - 1];
        }
        for ($j = $m; $j <= $high; ++$j) {
          $g = 0.0;
          for ($i = $m; $i <= $high; ++$i) {
            $g += $this->ort[$i] * $this->V[$i][$j];
          }

          // Double division avoids possible underflow
          $g = $g / $this->ort[$m] / $this->H[$m][$m - 1];
          for ($i = $m; $i <= $high; ++$i) {
            $this->V[$i][$j] += $g * $this->ort[$i];
          }
        }
      }
    }
  }

  /**
   *	Performs complex division.
   *
   *	@access private
   */
  private function cdiv($xr, $xi, $yr, $yi) {
    if (abs($yr) > abs($yi)) {
      $r = $yi / $yr;
      $d = $yr + $r * $yi;
      $this->cdivr = ($xr + $r * $xi) / $d;
      $this->cdivi = ($xi - $r * $xr) / $d;
    }
    else {
      $r = $yr / $yi;
      $d = $yi + $r * $yr;
      $this->cdivr = ($r * $xr + $xi) / $d;
      $this->cdivi = ($r * $xi - $xr) / $d;
    }
  }

  /**
   *	Nonsymmetric reduction from Hessenberg to real Schur form.
   *
   *	Code is derived from the Algol procedure hqr2,
   *	by Martin and Wilkinson, Handbook for Auto. Comp.,
   *	Vol.ii-Linear Algebra, and the corresponding
   *	Fortran subroutine in EISPACK.
   *
   *	@access private
   */
  private function hqr2() {

    //  Initialize
    $nn = $this->n;
    $n = $nn - 1;
    $low = 0;
    $high = $nn - 1;
    $eps = pow(2.0, -52.0);
    $exshift = 0.0;
    $p = $q = $r = $s = $z = 0;

    // Store roots isolated by balanc and compute matrix norm
    $norm = 0.0;
    for ($i = 0; $i < $nn; ++$i) {
      if ($i < $low or $i > $high) {
        $this->d[$i] = $this->H[$i][$i];
        $this->e[$i] = 0.0;
      }
      for ($j = max($i - 1, 0); $j < $nn; ++$j) {
        $norm = $norm + abs($this->H[$i][$j]);
      }
    }

    // Outer loop over eigenvalue index
    $iter = 0;
    while ($n >= $low) {

      // Look for single small sub-diagonal element
      $l = $n;
      while ($l > $low) {
        $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]);
        if ($s == 0.0) {
          $s = $norm;
        }
        if (abs($this->H[$l][$l - 1]) < $eps * $s) {
          break;
        }
        --$l;
      }

      // Check for convergence
      // One root found
      if ($l == $n) {
        $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
        $this->d[$n] = $this->H[$n][$n];
        $this->e[$n] = 0.0;
        --$n;
        $iter = 0;

        // Two roots found
      }
      else {
        if ($l == $n - 1) {
          $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
          $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0;
          $q = $p * $p + $w;
          $z = sqrt(abs($q));
          $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
          $this->H[$n - 1][$n - 1] = $this->H[$n - 1][$n - 1] + $exshift;
          $x = $this->H[$n][$n];

          // Real pair
          if ($q >= 0) {
            if ($p >= 0) {
              $z = $p + $z;
            }
            else {
              $z = $p - $z;
            }
            $this->d[$n - 1] = $x + $z;
            $this->d[$n] = $this->d[$n - 1];
            if ($z != 0.0) {
              $this->d[$n] = $x - $w / $z;
            }
            $this->e[$n - 1] = 0.0;
            $this->e[$n] = 0.0;
            $x = $this->H[$n][$n - 1];
            $s = abs($x) + abs($z);
            $p = $x / $s;
            $q = $z / $s;
            $r = sqrt($p * $p + $q * $q);
            $p = $p / $r;
            $q = $q / $r;

            // Row modification
            for ($j = $n - 1; $j < $nn; ++$j) {
              $z = $this->H[$n - 1][$j];
              $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j];
              $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
            }

            // Column modification
            for ($i = 0; $i <= n; ++$i) {
              $z = $this->H[$i][$n - 1];
              $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n];
              $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
            }

            // Accumulate transformations
            for ($i = $low; $i <= $high; ++$i) {
              $z = $this->V[$i][$n - 1];
              $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n];
              $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
            }

            // Complex pair
          }
          else {
            $this->d[$n - 1] = $x + $p;
            $this->d[$n] = $x + $p;
            $this->e[$n - 1] = $z;
            $this->e[$n] = -$z;
          }
          $n = $n - 2;
          $iter = 0;

          // No convergence yet
        }
        else {

          // Form shift
          $x = $this->H[$n][$n];
          $y = 0.0;
          $w = 0.0;
          if ($l < $n) {
            $y = $this->H[$n - 1][$n - 1];
            $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
          }

          // Wilkinson's original ad hoc shift
          if ($iter == 10) {
            $exshift += $x;
            for ($i = $low; $i <= $n; ++$i) {
              $this->H[$i][$i] -= $x;
            }
            $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]);
            $x = $y = 0.75 * $s;
            $w = -0.4375 * $s * $s;
          }

          // MATLAB's new ad hoc shift
          if ($iter == 30) {
            $s = ($y - $x) / 2.0;
            $s = $s * $s + $w;
            if ($s > 0) {
              $s = sqrt($s);
              if ($y < $x) {
                $s = -$s;
              }
              $s = $x - $w / (($y - $x) / 2.0 + $s);
              for ($i = $low; $i <= $n; ++$i) {
                $this->H[$i][$i] -= $s;
              }
              $exshift += $s;
              $x = $y = $w = 0.964;
            }
          }

          // Could check iteration count here.
          $iter = $iter + 1;

          // Look for two consecutive small sub-diagonal elements
          $m = $n - 2;
          while ($m >= $l) {
            $z = $this->H[$m][$m];
            $r = $x - $z;
            $s = $y - $z;
            $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1];
            $q = $this->H[$m + 1][$m + 1] - $z - $r - $s;
            $r = $this->H[$m + 2][$m + 1];
            $s = abs($p) + abs($q) + abs($r);
            $p = $p / $s;
            $q = $q / $s;
            $r = $r / $s;
            if ($m == $l) {
              break;
            }
            if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) < $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))) {
              break;
            }
            --$m;
          }
          for ($i = $m + 2; $i <= $n; ++$i) {
            $this->H[$i][$i - 2] = 0.0;
            if ($i > $m + 2) {
              $this->H[$i][$i - 3] = 0.0;
            }
          }

          // Double QR step involving rows l:n and columns m:n
          for ($k = $m; $k <= $n - 1; ++$k) {
            $notlast = $k != $n - 1;
            if ($k != $m) {
              $p = $this->H[$k][$k - 1];
              $q = $this->H[$k + 1][$k - 1];
              $r = $notlast ? $this->H[$k + 2][$k - 1] : 0.0;
              $x = abs($p) + abs($q) + abs($r);
              if ($x != 0.0) {
                $p = $p / $x;
                $q = $q / $x;
                $r = $r / $x;
              }
            }
            if ($x == 0.0) {
              break;
            }
            $s = sqrt($p * $p + $q * $q + $r * $r);
            if ($p < 0) {
              $s = -$s;
            }
            if ($s != 0) {
              if ($k != $m) {
                $this->H[$k][$k - 1] = -$s * $x;
              }
              elseif ($l != $m) {
                $this->H[$k][$k - 1] = -$this->H[$k][$k - 1];
              }
              $p = $p + $s;
              $x = $p / $s;
              $y = $q / $s;
              $z = $r / $s;
              $q = $q / $p;
              $r = $r / $p;

              // Row modification
              for ($j = $k; $j < $nn; ++$j) {
                $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j];
                if ($notlast) {
                  $p = $p + $r * $this->H[$k + 2][$j];
                  $this->H[$k + 2][$j] = $this->H[$k + 2][$j] - $p * $z;
                }
                $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
                $this->H[$k + 1][$j] = $this->H[$k + 1][$j] - $p * $y;
              }

              // Column modification
              for ($i = 0; $i <= min($n, $k + 3); ++$i) {
                $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1];
                if ($notlast) {
                  $p = $p + $z * $this->H[$i][$k + 2];
                  $this->H[$i][$k + 2] = $this->H[$i][$k + 2] - $p * $r;
                }
                $this->H[$i][$k] = $this->H[$i][$k] - $p;
                $this->H[$i][$k + 1] = $this->H[$i][$k + 1] - $p * $q;
              }

              // Accumulate transformations
              for ($i = $low; $i <= $high; ++$i) {
                $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1];
                if ($notlast) {
                  $p = $p + $z * $this->V[$i][$k + 2];
                  $this->V[$i][$k + 2] = $this->V[$i][$k + 2] - $p * $r;
                }
                $this->V[$i][$k] = $this->V[$i][$k] - $p;
                $this->V[$i][$k + 1] = $this->V[$i][$k + 1] - $p * $q;
              }
            }

            // ($s != 0)
          }

          // k loop
        }
      }

      // check convergence
    }

    // while ($n >= $low)
    // Backsubstitute to find vectors of upper triangular form
    if ($norm == 0.0) {
      return;
    }
    for ($n = $nn - 1; $n >= 0; --$n) {
      $p = $this->d[$n];
      $q = $this->e[$n];

      // Real vector
      if ($q == 0) {
        $l = $n;
        $this->H[$n][$n] = 1.0;
        for ($i = $n - 1; $i >= 0; --$i) {
          $w = $this->H[$i][$i] - $p;
          $r = 0.0;
          for ($j = $l; $j <= $n; ++$j) {
            $r = $r + $this->H[$i][$j] * $this->H[$j][$n];
          }
          if ($this->e[$i] < 0.0) {
            $z = $w;
            $s = $r;
          }
          else {
            $l = $i;
            if ($this->e[$i] == 0.0) {
              if ($w != 0.0) {
                $this->H[$i][$n] = -$r / $w;
              }
              else {
                $this->H[$i][$n] = -$r / ($eps * $norm);
              }

              // Solve real equations
            }
            else {
              $x = $this->H[$i][$i + 1];
              $y = $this->H[$i + 1][$i];
              $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
              $t = ($x * $s - $z * $r) / $q;
              $this->H[$i][$n] = $t;
              if (abs($x) > abs($z)) {
                $this->H[$i + 1][$n] = (-$r - $w * $t) / $x;
              }
              else {
                $this->H[$i + 1][$n] = (-$s - $y * $t) / $z;
              }
            }

            // Overflow control
            $t = abs($this->H[$i][$n]);
            if ($eps * $t * $t > 1) {
              for ($j = $i; $j <= $n; ++$j) {
                $this->H[$j][$n] = $this->H[$j][$n] / $t;
              }
            }
          }
        }

        // Complex vector
      }
      else {
        if ($q < 0) {
          $l = $n - 1;

          // Last vector component imaginary so matrix is triangular
          if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n])) {
            $this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1];
            $this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1];
          }
          else {
            $this
              ->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q);
            $this->H[$n - 1][$n - 1] = $this->cdivr;
            $this->H[$n - 1][$n] = $this->cdivi;
          }
          $this->H[$n][$n - 1] = 0.0;
          $this->H[$n][$n] = 1.0;
          for ($i = $n - 2; $i >= 0; --$i) {

            // double ra,sa,vr,vi;
            $ra = 0.0;
            $sa = 0.0;
            for ($j = $l; $j <= $n; ++$j) {
              $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n - 1];
              $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
            }
            $w = $this->H[$i][$i] - $p;
            if ($this->e[$i] < 0.0) {
              $z = $w;
              $r = $ra;
              $s = $sa;
            }
            else {
              $l = $i;
              if ($this->e[$i] == 0) {
                $this
                  ->cdiv(-$ra, -$sa, $w, $q);
                $this->H[$i][$n - 1] = $this->cdivr;
                $this->H[$i][$n] = $this->cdivi;
              }
              else {

                // Solve complex equations
                $x = $this->H[$i][$i + 1];
                $y = $this->H[$i + 1][$i];
                $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
                $vi = ($this->d[$i] - $p) * 2.0 * $q;
                if ($vr == 0.0 & $vi == 0.0) {
                  $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
                }
                $this
                  ->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
                $this->H[$i][$n - 1] = $this->cdivr;
                $this->H[$i][$n] = $this->cdivi;
                if (abs($x) > abs($z) + abs($q)) {
                  $this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x;
                  $this->H[$i + 1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x;
                }
                else {
                  $this
                    ->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q);
                  $this->H[$i + 1][$n - 1] = $this->cdivr;
                  $this->H[$i + 1][$n] = $this->cdivi;
                }
              }

              // Overflow control
              $t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n]));
              if ($eps * $t * $t > 1) {
                for ($j = $i; $j <= $n; ++$j) {
                  $this->H[$j][$n - 1] = $this->H[$j][$n - 1] / $t;
                  $this->H[$j][$n] = $this->H[$j][$n] / $t;
                }
              }
            }

            // end else
          }

          // end for
        }
      }

      // end else for complex case
    }

    // end for
    // Vectors of isolated roots
    for ($i = 0; $i < $nn; ++$i) {
      if ($i < $low | $i > $high) {
        for ($j = $i; $j < $nn; ++$j) {
          $this->V[$i][$j] = $this->H[$i][$j];
        }
      }
    }

    // Back transformation to get eigenvectors of original matrix
    for ($j = $nn - 1; $j >= $low; --$j) {
      for ($i = $low; $i <= $high; ++$i) {
        $z = 0.0;
        for ($k = $low; $k <= min($j, $high); ++$k) {
          $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
        }
        $this->V[$i][$j] = $z;
      }
    }
  }

  // end hqr2

  /**
   *	Constructor: Check for symmetry, then construct the eigenvalue decomposition
   *
   *	@access public
   *	@param A  Square matrix
   *	@return Structure to access D and V.
   */
  public function __construct($Arg) {
    $this->A = $Arg
      ->getArray();
    $this->n = $Arg
      ->getColumnDimension();
    $issymmetric = true;
    for ($j = 0; $j < $this->n & $issymmetric; ++$j) {
      for ($i = 0; $i < $this->n & $issymmetric; ++$i) {
        $issymmetric = $this->A[$i][$j] == $this->A[$j][$i];
      }
    }
    if ($issymmetric) {
      $this->V = $this->A;

      // Tridiagonalize.
      $this
        ->tred2();

      // Diagonalize.
      $this
        ->tql2();
    }
    else {
      $this->H = $this->A;
      $this->ort = array();

      // Reduce to Hessenberg form.
      $this
        ->orthes();

      // Reduce Hessenberg to real Schur form.
      $this
        ->hqr2();
    }
  }

  /**
   *	Return the eigenvector matrix
   *
   *	@access public
   *	@return V
   */
  public function getV() {
    return new Matrix($this->V, $this->n, $this->n);
  }

  /**
   *	Return the real parts of the eigenvalues
   *
   *	@access public
   *	@return real(diag(D))
   */
  public function getRealEigenvalues() {
    return $this->d;
  }

  /**
   *	Return the imaginary parts of the eigenvalues
   *
   *	@access public
   *	@return imag(diag(D))
   */
  public function getImagEigenvalues() {
    return $this->e;
  }

  /**
   *	Return the block diagonal eigenvalue matrix
   *
   *	@access public
   *	@return D
   */
  public function getD() {
    for ($i = 0; $i < $this->n; ++$i) {
      $D[$i] = array_fill(0, $this->n, 0.0);
      $D[$i][$i] = $this->d[$i];
      if ($this->e[$i] == 0) {
        continue;
      }
      $o = $this->e[$i] > 0 ? $i + 1 : $i - 1;
      $D[$i][$o] = $this->e[$i];
    }
    return new Matrix($D);
  }

}

//	class EigenvalueDecomposition

Classes

Name	Description
EigenvalueDecomposition	@package JAMA