private function EigenvalueDecomposition::hqr2 in Loft Data Grids 7.2

Same name and namespace in other branches

1. 6.2 vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php \EigenvalueDecomposition::hqr2()

* 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

1 call to EigenvalueDecomposition::hqr2()

EigenvalueDecomposition::__construct in vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php

* Constructor: Check for symmetry, then construct the eigenvalue decomposition * * @access public *

File

vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php, line 398

Class

EigenvalueDecomposition

@package JAMA

Code

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;
    }
  }
}