Search moodle.org's
Developer Documentation

See Release Notes

  • Bug fixes for general core bugs in 3.11.x will end 14 Nov 2022 (12 months plus 6 months extension).
  • Bug fixes for security issues in 3.11.x will end 13 Nov 2023 (18 months plus 12 months extension).
  • PHP version: minimum PHP 7.3.0 Note: minimum PHP version has increased since Moodle 3.10. PHP 7.4.x is supported too.

Differences Between: [Versions 311 and 400] [Versions 311 and 401] [Versions 311 and 402] [Versions 311 and 403]

   1  <?php
   2  
   3  declare(strict_types=1);
   4  
   5  /**
   6   * Class to obtain eigenvalues and eigenvectors of a real matrix.
   7   *
   8   * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D
   9   * is diagonal and the eigenvector matrix V is orthogonal (i.e.
  10   * A = V.times(D.times(V.transpose())) and V.times(V.transpose())
  11   * equals the identity matrix).
  12   *
  13   * If A is not symmetric, then the eigenvalue matrix D is block diagonal
  14   * with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
  15   * lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda].  The
  16   * columns of V represent the eigenvectors in the sense that A*V = V*D,
  17   * i.e. A.times(V) equals V.times(D).  The matrix V may be badly
  18   * conditioned, or even singular, so the validity of the equation
  19   * A = V*D*inverse(V) depends upon V.cond().
  20   *
  21   * @author Paul Meagher
  22   * @license PHP v3.0
  23   *
  24   * @version 1.1
  25   *
  26   *  Slightly changed to adapt the original code to PHP-ML library
  27   *  @date 2017/04/11
  28   *
  29   *  @author Mustafa Karabulut
  30   */
  31  
  32  namespace Phpml\Math\LinearAlgebra;
  33  
  34  use Phpml\Math\Matrix;
  35  
  36  class EigenvalueDecomposition
  37  {
  38      /**
  39       * Row and column dimension (square matrix).
  40       *
  41       * @var int
  42       */
  43      private $n;
  44  
  45      /**
  46       * Arrays for internal storage of eigenvalues.
  47       *
  48       * @var array
  49       */
  50      private $d = [];
  51  
  52      /**
  53       * @var array
  54       */
  55      private $e = [];
  56  
  57      /**
  58       * Array for internal storage of eigenvectors.
  59       *
  60       * @var array
  61       */
  62      private $V = [];
  63  
  64      /**
  65       * Array for internal storage of nonsymmetric Hessenberg form.
  66       *
  67       * @var array
  68       */
  69      private $H = [];
  70  
  71      /**
  72       * Working storage for nonsymmetric algorithm.
  73       *
  74       * @var array
  75       */
  76      private $ort = [];
  77  
  78      /**
  79       * Used for complex scalar division.
  80       *
  81       * @var float
  82       */
  83      private $cdivr;
  84  
  85      /**
  86       * @var float
  87       */
  88      private $cdivi;
  89  
  90      /**
  91       * Constructor: Check for symmetry, then construct the eigenvalue decomposition
  92       */
  93      public function __construct(array $arg)
  94      {
  95          $this->n = count($arg[0]);
  96          $symmetric = true;
  97  
  98          for ($j = 0; ($j < $this->n) & $symmetric; ++$j) {
  99              for ($i = 0; ($i < $this->n) & $symmetric; ++$i) {
 100                  $symmetric = $arg[$i][$j] == $arg[$j][$i];
 101              }
 102          }
 103  
 104          if ($symmetric) {
 105              $this->V = $arg;
 106              // Tridiagonalize.
 107              $this->tred2();
 108              // Diagonalize.
 109              $this->tql2();
 110          } else {
 111              $this->H = $arg;
 112              $this->ort = [];
 113              // Reduce to Hessenberg form.
 114              $this->orthes();
 115              // Reduce Hessenberg to real Schur form.
 116              $this->hqr2();
 117          }
 118      }
 119  
 120      /**
 121       * Return the eigenvector matrix
 122       */
 123      public function getEigenvectors(): array
 124      {
 125          $vectors = $this->V;
 126  
 127          // Always return the eigenvectors of length 1.0
 128          $vectors = new Matrix($vectors);
 129          $vectors = array_map(function ($vect) {
 130              $sum = 0;
 131              $count = count($vect);
 132              for ($i = 0; $i < $count; ++$i) {
 133                  $sum += $vect[$i] ** 2;
 134              }
 135  
 136              $sum **= .5;
 137              for ($i = 0; $i < $count; ++$i) {
 138                  $vect[$i] /= $sum;
 139              }
 140  
 141              return $vect;
 142          }, $vectors->transpose()->toArray());
 143  
 144          return $vectors;
 145      }
 146  
 147      /**
 148       * Return the real parts of the eigenvalues<br>
 149       *  d = real(diag(D));
 150       */
 151      public function getRealEigenvalues(): array
 152      {
 153          return $this->d;
 154      }
 155  
 156      /**
 157       * Return the imaginary parts of the eigenvalues <br>
 158       *  d = imag(diag(D))
 159       */
 160      public function getImagEigenvalues(): array
 161      {
 162          return $this->e;
 163      }
 164  
 165      /**
 166       * Return the block diagonal eigenvalue matrix
 167       */
 168      public function getDiagonalEigenvalues(): array
 169      {
 170          $D = [];
 171  
 172          for ($i = 0; $i < $this->n; ++$i) {
 173              $D[$i] = array_fill(0, $this->n, 0.0);
 174              $D[$i][$i] = $this->d[$i];
 175              if ($this->e[$i] == 0) {
 176                  continue;
 177              }
 178  
 179              $o = $this->e[$i] > 0 ? $i + 1 : $i - 1;
 180              $D[$i][$o] = $this->e[$i];
 181          }
 182  
 183          return $D;
 184      }
 185  
 186      /**
 187       * Symmetric Householder reduction to tridiagonal form.
 188       */
 189      private function tred2(): void
 190      {
 191          //  This is derived from the Algol procedures tred2 by
 192          //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
 193          //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
 194          //  Fortran subroutine in EISPACK.
 195          $this->d = $this->V[$this->n - 1];
 196          // Householder reduction to tridiagonal form.
 197          for ($i = $this->n - 1; $i > 0; --$i) {
 198              $i_ = $i - 1;
 199              // Scale to avoid under/overflow.
 200              $h = $scale = 0.0;
 201              $scale += array_sum(array_map('abs', $this->d));
 202              if ($scale == 0.0) {
 203                  $this->e[$i] = $this->d[$i_];
 204                  $this->d = array_slice($this->V[$i_], 0, $this->n - 1);
 205                  for ($j = 0; $j < $i; ++$j) {
 206                      $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
 207                  }
 208              } else {
 209                  // Generate Householder vector.
 210                  for ($k = 0; $k < $i; ++$k) {
 211                      $this->d[$k] /= $scale;
 212                      $h += $this->d[$k] ** 2;
 213                  }
 214  
 215                  $f = $this->d[$i_];
 216                  $g = $h ** .5;
 217                  if ($f > 0) {
 218                      $g = -$g;
 219                  }
 220  
 221                  $this->e[$i] = $scale * $g;
 222                  $h -= $f * $g;
 223                  $this->d[$i_] = $f - $g;
 224  
 225                  for ($j = 0; $j < $i; ++$j) {
 226                      $this->e[$j] = 0.0;
 227                  }
 228  
 229                  // Apply similarity transformation to remaining columns.
 230                  for ($j = 0; $j < $i; ++$j) {
 231                      $f = $this->d[$j];
 232                      $this->V[$j][$i] = $f;
 233                      $g = $this->e[$j] + $this->V[$j][$j] * $f;
 234  
 235                      for ($k = $j + 1; $k <= $i_; ++$k) {
 236                          $g += $this->V[$k][$j] * $this->d[$k];
 237                          $this->e[$k] += $this->V[$k][$j] * $f;
 238                      }
 239  
 240                      $this->e[$j] = $g;
 241                  }
 242  
 243                  $f = 0.0;
 244  
 245                  if ($h == 0.0) {
 246                      $h = 1e-32;
 247                  }
 248  
 249                  for ($j = 0; $j < $i; ++$j) {
 250                      $this->e[$j] /= $h;
 251                      $f += $this->e[$j] * $this->d[$j];
 252                  }
 253  
 254                  $hh = $f / (2 * $h);
 255                  for ($j = 0; $j < $i; ++$j) {
 256                      $this->e[$j] -= $hh * $this->d[$j];
 257                  }
 258  
 259                  for ($j = 0; $j < $i; ++$j) {
 260                      $f = $this->d[$j];
 261                      $g = $this->e[$j];
 262                      for ($k = $j; $k <= $i_; ++$k) {
 263                          $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
 264                      }
 265  
 266                      $this->d[$j] = $this->V[$i - 1][$j];
 267                      $this->V[$i][$j] = 0.0;
 268                  }
 269              }
 270  
 271              $this->d[$i] = $h;
 272          }
 273  
 274          // Accumulate transformations.
 275          for ($i = 0; $i < $this->n - 1; ++$i) {
 276              $this->V[$this->n - 1][$i] = $this->V[$i][$i];
 277              $this->V[$i][$i] = 1.0;
 278              $h = $this->d[$i + 1];
 279              if ($h != 0.0) {
 280                  for ($k = 0; $k <= $i; ++$k) {
 281                      $this->d[$k] = $this->V[$k][$i + 1] / $h;
 282                  }
 283  
 284                  for ($j = 0; $j <= $i; ++$j) {
 285                      $g = 0.0;
 286                      for ($k = 0; $k <= $i; ++$k) {
 287                          $g += $this->V[$k][$i + 1] * $this->V[$k][$j];
 288                      }
 289  
 290                      for ($k = 0; $k <= $i; ++$k) {
 291                          $this->V[$k][$j] -= $g * $this->d[$k];
 292                      }
 293                  }
 294              }
 295  
 296              for ($k = 0; $k <= $i; ++$k) {
 297                  $this->V[$k][$i + 1] = 0.0;
 298              }
 299          }
 300  
 301          $this->d = $this->V[$this->n - 1];
 302          $this->V[$this->n - 1] = array_fill(0, $this->n, 0.0);
 303          $this->V[$this->n - 1][$this->n - 1] = 1.0;
 304          $this->e[0] = 0.0;
 305      }
 306  
 307      /**
 308       * Symmetric tridiagonal QL algorithm.
 309       *
 310       * This is derived from the Algol procedures tql2, by
 311       * Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
 312       * Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
 313       * Fortran subroutine in EISPACK.
 314       */
 315      private function tql2(): void
 316      {
 317          for ($i = 1; $i < $this->n; ++$i) {
 318              $this->e[$i - 1] = $this->e[$i];
 319          }
 320  
 321          $this->e[$this->n - 1] = 0.0;
 322          $f = 0.0;
 323          $tst1 = 0.0;
 324          $eps = 2.0 ** -52.0;
 325  
 326          for ($l = 0; $l < $this->n; ++$l) {
 327              // Find small subdiagonal element
 328              $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
 329              $m = $l;
 330              while ($m < $this->n) {
 331                  if (abs($this->e[$m]) <= $eps * $tst1) {
 332                      break;
 333                  }
 334  
 335                  ++$m;
 336              }
 337  
 338              // If m == l, $this->d[l] is an eigenvalue,
 339              // otherwise, iterate.
 340              if ($m > $l) {
 341                  do {
 342                      // Compute implicit shift
 343                      $g = $this->d[$l];
 344                      $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]);
 345                      $r = hypot($p, 1.0);
 346                      if ($p < 0) {
 347                          $r *= -1;
 348                      }
 349  
 350                      $this->d[$l] = $this->e[$l] / ($p + $r);
 351                      $this->d[$l + 1] = $this->e[$l] * ($p + $r);
 352                      $dl1 = $this->d[$l + 1];
 353                      $h = $g - $this->d[$l];
 354                      for ($i = $l + 2; $i < $this->n; ++$i) {
 355                          $this->d[$i] -= $h;
 356                      }
 357  
 358                      $f += $h;
 359                      // Implicit QL transformation.
 360                      $p = $this->d[$m];
 361                      $c = 1.0;
 362                      $c2 = $c3 = $c;
 363                      $el1 = $this->e[$l + 1];
 364                      $s = $s2 = 0.0;
 365                      for ($i = $m - 1; $i >= $l; --$i) {
 366                          $c3 = $c2;
 367                          $c2 = $c;
 368                          $s2 = $s;
 369                          $g = $c * $this->e[$i];
 370                          $h = $c * $p;
 371                          $r = hypot($p, $this->e[$i]);
 372                          $this->e[$i + 1] = $s * $r;
 373                          $s = $this->e[$i] / $r;
 374                          $c = $p / $r;
 375                          $p = $c * $this->d[$i] - $s * $g;
 376                          $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]);
 377                          // Accumulate transformation.
 378                          for ($k = 0; $k < $this->n; ++$k) {
 379                              $h = $this->V[$k][$i + 1];
 380                              $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h;
 381                              $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
 382                          }
 383                      }
 384  
 385                      $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
 386                      $this->e[$l] = $s * $p;
 387                      $this->d[$l] = $c * $p;
 388                      // Check for convergence.
 389                  } while (abs($this->e[$l]) > $eps * $tst1);
 390              }
 391  
 392              $this->d[$l] += $f;
 393              $this->e[$l] = 0.0;
 394          }
 395  
 396          // Sort eigenvalues and corresponding vectors.
 397          for ($i = 0; $i < $this->n - 1; ++$i) {
 398              $k = $i;
 399              $p = $this->d[$i];
 400              for ($j = $i + 1; $j < $this->n; ++$j) {
 401                  if ($this->d[$j] < $p) {
 402                      $k = $j;
 403                      $p = $this->d[$j];
 404                  }
 405              }
 406  
 407              if ($k != $i) {
 408                  $this->d[$k] = $this->d[$i];
 409                  $this->d[$i] = $p;
 410                  for ($j = 0; $j < $this->n; ++$j) {
 411                      $p = $this->V[$j][$i];
 412                      $this->V[$j][$i] = $this->V[$j][$k];
 413                      $this->V[$j][$k] = $p;
 414                  }
 415              }
 416          }
 417      }
 418  
 419      /**
 420       * Nonsymmetric reduction to Hessenberg form.
 421       *
 422       * This is derived from the Algol procedures orthes and ortran,
 423       * by Martin and Wilkinson, Handbook for Auto. Comp.,
 424       * Vol.ii-Linear Algebra, and the corresponding
 425       * Fortran subroutines in EISPACK.
 426       */
 427      private function orthes(): void
 428      {
 429          $low = 0;
 430          $high = $this->n - 1;
 431  
 432          for ($m = $low + 1; $m <= $high - 1; ++$m) {
 433              // Scale column.
 434              $scale = 0.0;
 435              for ($i = $m; $i <= $high; ++$i) {
 436                  $scale += abs($this->H[$i][$m - 1]);
 437              }
 438  
 439              if ($scale != 0.0) {
 440                  // Compute Householder transformation.
 441                  $h = 0.0;
 442                  for ($i = $high; $i >= $m; --$i) {
 443                      $this->ort[$i] = $this->H[$i][$m - 1] / $scale;
 444                      $h += $this->ort[$i] * $this->ort[$i];
 445                  }
 446  
 447                  $g = $h ** .5;
 448                  if ($this->ort[$m] > 0) {
 449                      $g *= -1;
 450                  }
 451  
 452                  $h -= $this->ort[$m] * $g;
 453                  $this->ort[$m] -= $g;
 454                  // Apply Householder similarity transformation
 455                  // H = (I -u * u' / h) * H * (I -u * u') / h)
 456                  for ($j = $m; $j < $this->n; ++$j) {
 457                      $f = 0.0;
 458                      for ($i = $high; $i >= $m; --$i) {
 459                          $f += $this->ort[$i] * $this->H[$i][$j];
 460                      }
 461  
 462                      $f /= $h;
 463                      for ($i = $m; $i <= $high; ++$i) {
 464                          $this->H[$i][$j] -= $f * $this->ort[$i];
 465                      }
 466                  }
 467  
 468                  for ($i = 0; $i <= $high; ++$i) {
 469                      $f = 0.0;
 470                      for ($j = $high; $j >= $m; --$j) {
 471                          $f += $this->ort[$j] * $this->H[$i][$j];
 472                      }
 473  
 474                      $f /= $h;
 475                      for ($j = $m; $j <= $high; ++$j) {
 476                          $this->H[$i][$j] -= $f * $this->ort[$j];
 477                      }
 478                  }
 479  
 480                  $this->ort[$m] = $scale * $this->ort[$m];
 481                  $this->H[$m][$m - 1] = $scale * $g;
 482              }
 483          }
 484  
 485          // Accumulate transformations (Algol's ortran).
 486          for ($i = 0; $i < $this->n; ++$i) {
 487              for ($j = 0; $j < $this->n; ++$j) {
 488                  $this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
 489              }
 490          }
 491  
 492          for ($m = $high - 1; $m >= $low + 1; --$m) {
 493              if ($this->H[$m][$m - 1] != 0.0) {
 494                  for ($i = $m + 1; $i <= $high; ++$i) {
 495                      $this->ort[$i] = $this->H[$i][$m - 1];
 496                  }
 497  
 498                  for ($j = $m; $j <= $high; ++$j) {
 499                      $g = 0.0;
 500                      for ($i = $m; $i <= $high; ++$i) {
 501                          $g += $this->ort[$i] * $this->V[$i][$j];
 502                      }
 503  
 504                      // Double division avoids possible underflow
 505                      $g = ($g / $this->ort[$m]) / $this->H[$m][$m - 1];
 506                      for ($i = $m; $i <= $high; ++$i) {
 507                          $this->V[$i][$j] += $g * $this->ort[$i];
 508                      }
 509                  }
 510              }
 511          }
 512      }
 513  
 514      /**
 515       * Performs complex division.
 516       *
 517       * @param int|float $xr
 518       * @param int|float $xi
 519       * @param int|float $yr
 520       * @param int|float $yi
 521       */
 522      private function cdiv($xr, $xi, $yr, $yi): void
 523      {
 524          if (abs($yr) > abs($yi)) {
 525              $r = $yi / $yr;
 526              $d = $yr + $r * $yi;
 527              $this->cdivr = ($xr + $r * $xi) / $d;
 528              $this->cdivi = ($xi - $r * $xr) / $d;
 529          } else {
 530              $r = $yr / $yi;
 531              $d = $yi + $r * $yr;
 532              $this->cdivr = ($r * $xr + $xi) / $d;
 533              $this->cdivi = ($r * $xi - $xr) / $d;
 534          }
 535      }
 536  
 537      /**
 538       * Nonsymmetric reduction from Hessenberg to real Schur form.
 539       *
 540       * Code is derived from the Algol procedure hqr2,
 541       * by Martin and Wilkinson, Handbook for Auto. Comp.,
 542       * Vol.ii-Linear Algebra, and the corresponding
 543       * Fortran subroutine in EISPACK.
 544       */
 545      private function hqr2(): void
 546      {
 547          //  Initialize
 548          $nn = $this->n;
 549          $n = $nn - 1;
 550          $low = 0;
 551          $high = $nn - 1;
 552          $eps = 2.0 ** -52.0;
 553          $exshift = 0.0;
 554          $p = $q = $r = $s = $z = 0;
 555          // Store roots isolated by balanc and compute matrix norm
 556          $norm = 0.0;
 557  
 558          for ($i = 0; $i < $nn; ++$i) {
 559              if (($i < $low) or ($i > $high)) {
 560                  $this->d[$i] = $this->H[$i][$i];
 561                  $this->e[$i] = 0.0;
 562              }
 563  
 564              for ($j = max($i - 1, 0); $j < $nn; ++$j) {
 565                  $norm += abs($this->H[$i][$j]);
 566              }
 567          }
 568  
 569          // Outer loop over eigenvalue index
 570          $iter = 0;
 571          while ($n >= $low) {
 572              // Look for single small sub-diagonal element
 573              $l = $n;
 574              while ($l > $low) {
 575                  $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]);
 576                  if ($s == 0.0) {
 577                      $s = $norm;
 578                  }
 579  
 580                  if (abs($this->H[$l][$l - 1]) < $eps * $s) {
 581                      break;
 582                  }
 583  
 584                  --$l;
 585              }
 586  
 587              // Check for convergence
 588              // One root found
 589              if ($l == $n) {
 590                  $this->H[$n][$n] += $exshift;
 591                  $this->d[$n] = $this->H[$n][$n];
 592                  $this->e[$n] = 0.0;
 593                  --$n;
 594                  $iter = 0;
 595              // Two roots found
 596              } elseif ($l == $n - 1) {
 597                  $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
 598                  $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0;
 599                  $q = $p * $p + $w;
 600                  $z = abs($q) ** .5;
 601                  $this->H[$n][$n] += $exshift;
 602                  $this->H[$n - 1][$n - 1] += $exshift;
 603                  $x = $this->H[$n][$n];
 604                  // Real pair
 605                  if ($q >= 0) {
 606                      if ($p >= 0) {
 607                          $z = $p + $z;
 608                      } else {
 609                          $z = $p - $z;
 610                      }
 611  
 612                      $this->d[$n - 1] = $x + $z;
 613                      $this->d[$n] = $this->d[$n - 1];
 614                      if ($z != 0.0) {
 615                          $this->d[$n] = $x - $w / $z;
 616                      }
 617  
 618                      $this->e[$n - 1] = 0.0;
 619                      $this->e[$n] = 0.0;
 620                      $x = $this->H[$n][$n - 1];
 621                      $s = abs($x) + abs($z);
 622                      $p = $x / $s;
 623                      $q = $z / $s;
 624                      $r = ($p * $p + $q * $q) ** .5;
 625                      $p /= $r;
 626                      $q /= $r;
 627                      // Row modification
 628                      for ($j = $n - 1; $j < $nn; ++$j) {
 629                          $z = $this->H[$n - 1][$j];
 630                          $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j];
 631                          $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
 632                      }
 633  
 634                      // Column modification
 635                      for ($i = 0; $i <= $n; ++$i) {
 636                          $z = $this->H[$i][$n - 1];
 637                          $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n];
 638                          $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
 639                      }
 640  
 641                      // Accumulate transformations
 642                      for ($i = $low; $i <= $high; ++$i) {
 643                          $z = $this->V[$i][$n - 1];
 644                          $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n];
 645                          $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
 646                      }
 647  
 648                      // Complex pair
 649                  } else {
 650                      $this->d[$n - 1] = $x + $p;
 651                      $this->d[$n] = $x + $p;
 652                      $this->e[$n - 1] = $z;
 653                      $this->e[$n] = -$z;
 654                  }
 655  
 656                  $n -= 2;
 657                  $iter = 0;
 658              // No convergence yet
 659              } else {
 660                  // Form shift
 661                  $x = $this->H[$n][$n];
 662                  $y = 0.0;
 663                  $w = 0.0;
 664                  if ($l < $n) {
 665                      $y = $this->H[$n - 1][$n - 1];
 666                      $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
 667                  }
 668  
 669                  // Wilkinson's original ad hoc shift
 670                  if ($iter == 10) {
 671                      $exshift += $x;
 672                      for ($i = $low; $i <= $n; ++$i) {
 673                          $this->H[$i][$i] -= $x;
 674                      }
 675  
 676                      $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]);
 677                      $x = $y = 0.75 * $s;
 678                      $w = -0.4375 * $s * $s;
 679                  }
 680  
 681                  // MATLAB's new ad hoc shift
 682                  if ($iter == 30) {
 683                      $s = ($y - $x) / 2.0;
 684                      $s *= $s + $w;
 685                      if ($s > 0) {
 686                          $s **= .5;
 687                          if ($y < $x) {
 688                              $s = -$s;
 689                          }
 690  
 691                          $s = $x - $w / (($y - $x) / 2.0 + $s);
 692                          for ($i = $low; $i <= $n; ++$i) {
 693                              $this->H[$i][$i] -= $s;
 694                          }
 695  
 696                          $exshift += $s;
 697                          $x = $y = $w = 0.964;
 698                      }
 699                  }
 700  
 701                  // Could check iteration count here.
 702                  ++$iter;
 703                  // Look for two consecutive small sub-diagonal elements
 704                  $m = $n - 2;
 705                  while ($m >= $l) {
 706                      $z = $this->H[$m][$m];
 707                      $r = $x - $z;
 708                      $s = $y - $z;
 709                      $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1];
 710                      $q = $this->H[$m + 1][$m + 1] - $z - $r - $s;
 711                      $r = $this->H[$m + 2][$m + 1];
 712                      $s = abs($p) + abs($q) + abs($r);
 713                      $p /= $s;
 714                      $q /= $s;
 715                      $r /= $s;
 716                      if ($m == $l) {
 717                          break;
 718                      }
 719  
 720                      if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) <
 721                          $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))) {
 722                          break;
 723                      }
 724  
 725                      --$m;
 726                  }
 727  
 728                  for ($i = $m + 2; $i <= $n; ++$i) {
 729                      $this->H[$i][$i - 2] = 0.0;
 730                      if ($i > $m + 2) {
 731                          $this->H[$i][$i - 3] = 0.0;
 732                      }
 733                  }
 734  
 735                  // Double QR step involving rows l:n and columns m:n
 736                  for ($k = $m; $k <= $n - 1; ++$k) {
 737                      $notlast = ($k != $n - 1);
 738                      if ($k != $m) {
 739                          $p = $this->H[$k][$k - 1];
 740                          $q = $this->H[$k + 1][$k - 1];
 741                          $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0);
 742                          $x = abs($p) + abs($q) + abs($r);
 743                          if ($x != 0.0) {
 744                              $p /= $x;
 745                              $q /= $x;
 746                              $r /= $x;
 747                          }
 748                      }
 749  
 750                      if ($x == 0.0) {
 751                          break;
 752                      }
 753  
 754                      $s = ($p * $p + $q * $q + $r * $r) ** .5;
 755                      if ($p < 0) {
 756                          $s = -$s;
 757                      }
 758  
 759                      if ($s != 0) {
 760                          if ($k != $m) {
 761                              $this->H[$k][$k - 1] = -$s * $x;
 762                          } elseif ($l != $m) {
 763                              $this->H[$k][$k - 1] = -$this->H[$k][$k - 1];
 764                          }
 765  
 766                          $p += $s;
 767                          $x = $p / $s;
 768                          $y = $q / $s;
 769                          $z = $r / $s;
 770                          $q /= $p;
 771                          $r /= $p;
 772                          // Row modification
 773                          for ($j = $k; $j < $nn; ++$j) {
 774                              $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j];
 775                              if ($notlast) {
 776                                  $p += $r * $this->H[$k + 2][$j];
 777                                  $this->H[$k + 2][$j] -= $p * $z;
 778                              }
 779  
 780                              $this->H[$k][$j] -= $p * $x;
 781                              $this->H[$k + 1][$j] -= $p * $y;
 782                          }
 783  
 784                          // Column modification
 785                          for ($i = 0; $i <= min($n, $k + 3); ++$i) {
 786                              $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1];
 787                              if ($notlast) {
 788                                  $p += $z * $this->H[$i][$k + 2];
 789                                  $this->H[$i][$k + 2] -= $p * $r;
 790                              }
 791  
 792                              $this->H[$i][$k] -= $p;
 793                              $this->H[$i][$k + 1] -= $p * $q;
 794                          }
 795  
 796                          // Accumulate transformations
 797                          for ($i = $low; $i <= $high; ++$i) {
 798                              $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1];
 799                              if ($notlast) {
 800                                  $p += $z * $this->V[$i][$k + 2];
 801                                  $this->V[$i][$k + 2] -= $p * $r;
 802                              }
 803  
 804                              $this->V[$i][$k] -= $p;
 805                              $this->V[$i][$k + 1] -= $p * $q;
 806                          }
 807                      }  // ($s != 0)
 808                  }  // k loop
 809              }  // check convergence
 810          }  // while ($n >= $low)
 811  
 812          // Backsubstitute to find vectors of upper triangular form
 813          if ($norm == 0.0) {
 814              return;
 815          }
 816  
 817          for ($n = $nn - 1; $n >= 0; --$n) {
 818              $p = $this->d[$n];
 819              $q = $this->e[$n];
 820              // Real vector
 821              if ($q == 0) {
 822                  $l = $n;
 823                  $this->H[$n][$n] = 1.0;
 824                  for ($i = $n - 1; $i >= 0; --$i) {
 825                      $w = $this->H[$i][$i] - $p;
 826                      $r = 0.0;
 827                      for ($j = $l; $j <= $n; ++$j) {
 828                          $r += $this->H[$i][$j] * $this->H[$j][$n];
 829                      }
 830  
 831                      if ($this->e[$i] < 0.0) {
 832                          $z = $w;
 833                          $s = $r;
 834                      } else {
 835                          $l = $i;
 836                          if ($this->e[$i] == 0.0) {
 837                              if ($w != 0.0) {
 838                                  $this->H[$i][$n] = -$r / $w;
 839                              } else {
 840                                  $this->H[$i][$n] = -$r / ($eps * $norm);
 841                              }
 842  
 843                              // Solve real equations
 844                          } else {
 845                              $x = $this->H[$i][$i + 1];
 846                              $y = $this->H[$i + 1][$i];
 847                              $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
 848                              $t = ($x * $s - $z * $r) / $q;
 849                              $this->H[$i][$n] = $t;
 850                              if (abs($x) > abs($z)) {
 851                                  $this->H[$i + 1][$n] = (-$r - $w * $t) / $x;
 852                              } else {
 853                                  $this->H[$i + 1][$n] = (-$s - $y * $t) / $z;
 854                              }
 855                          }
 856  
 857                          // Overflow control
 858                          $t = abs($this->H[$i][$n]);
 859                          if (($eps * $t) * $t > 1) {
 860                              for ($j = $i; $j <= $n; ++$j) {
 861                                  $this->H[$j][$n] /= $t;
 862                              }
 863                          }
 864                      }
 865                  }
 866  
 867                  // Complex vector
 868              } elseif ($q < 0) {
 869                  $l = $n - 1;
 870                  // Last vector component imaginary so matrix is triangular
 871                  if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n])) {
 872                      $this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1];
 873                      $this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1];
 874                  } else {
 875                      $this->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q);
 876                      $this->H[$n - 1][$n - 1] = $this->cdivr;
 877                      $this->H[$n - 1][$n] = $this->cdivi;
 878                  }
 879  
 880                  $this->H[$n][$n - 1] = 0.0;
 881                  $this->H[$n][$n] = 1.0;
 882                  for ($i = $n - 2; $i >= 0; --$i) {
 883                      // double ra,sa,vr,vi;
 884                      $ra = 0.0;
 885                      $sa = 0.0;
 886                      for ($j = $l; $j <= $n; ++$j) {
 887                          $ra += $this->H[$i][$j] * $this->H[$j][$n - 1];
 888                          $sa += $this->H[$i][$j] * $this->H[$j][$n];
 889                      }
 890  
 891                      $w = $this->H[$i][$i] - $p;
 892                      if ($this->e[$i] < 0.0) {
 893                          $z = $w;
 894                          $r = $ra;
 895                          $s = $sa;
 896                      } else {
 897                          $l = $i;
 898                          if ($this->e[$i] == 0) {
 899                              $this->cdiv(-$ra, -$sa, $w, $q);
 900                              $this->H[$i][$n - 1] = $this->cdivr;
 901                              $this->H[$i][$n] = $this->cdivi;
 902                          } else {
 903                              // Solve complex equations
 904                              $x = $this->H[$i][$i + 1];
 905                              $y = $this->H[$i + 1][$i];
 906                              $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
 907                              $vi = ($this->d[$i] - $p) * 2.0 * $q;
 908                              if ($vr == 0.0 && $vi == 0.0) {
 909                                  $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
 910                              }
 911  
 912                              $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
 913                              $this->H[$i][$n - 1] = $this->cdivr;
 914                              $this->H[$i][$n] = $this->cdivi;
 915                              if (abs($x) > (abs($z) + abs($q))) {
 916                                  $this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x;
 917                                  $this->H[$i + 1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x;
 918                              } else {
 919                                  $this->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q);
 920                                  $this->H[$i + 1][$n - 1] = $this->cdivr;
 921                                  $this->H[$i + 1][$n] = $this->cdivi;
 922                              }
 923                          }
 924  
 925                          // Overflow control
 926                          $t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n]));
 927                          if (($eps * $t) * $t > 1) {
 928                              for ($j = $i; $j <= $n; ++$j) {
 929                                  $this->H[$j][$n - 1] /= $t;
 930                                  $this->H[$j][$n] /= $t;
 931                              }
 932                          }
 933                      } // end else
 934                  } // end for
 935              } // end else for complex case
 936          } // end for
 937  
 938          // Vectors of isolated roots
 939          for ($i = 0; $i < $nn; ++$i) {
 940              if ($i < $low || $i > $high) {
 941                  for ($j = $i; $j < $nn; ++$j) {
 942                      $this->V[$i][$j] = $this->H[$i][$j];
 943                  }
 944              }
 945          }
 946  
 947          // Back transformation to get eigenvectors of original matrix
 948          for ($j = $nn - 1; $j >= $low; --$j) {
 949              for ($i = $low; $i <= $high; ++$i) {
 950                  $z = 0.0;
 951                  for ($k = $low; $k <= min($j, $high); ++$k) {
 952                      $z += $this->V[$i][$k] * $this->H[$k][$j];
 953                  }
 954  
 955                  $this->V[$i][$j] = $z;
 956              }
 957          }
 958      }
 959  }