81                $this->d = $this->V[$this->n-1];
 
   83                for ($i = $this->n-1; $i > 0; --$i) {
 
   87                        $scale += array_sum(array_map(abs, $this->d));
 
   89                                $this->
e[$i] = $this->d[$i_];
 
   90                                $this->d = array_slice($this->V[$i_], 0, $i_);
 
   91                                for ($j = 0; $j < $i; ++$j) {
 
   92                                        $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
 
   96                                for ($k = 0; $k < $i; ++$k) {
 
   97                                        $this->d[$k] /= $scale;
 
   98                                        $h += pow($this->d[$k], 2);
 
  105                                $this->
e[$i] = $scale * $g;
 
  107                                $this->d[$i_] = $f - $g;
 
  108                                for ($j = 0; $j < $i; ++$j) {
 
  112                                for ($j = 0; $j < $i; ++$j) {
 
  114                                        $this->V[$j][$i] = $f;
 
  115                                        $g = $this->
e[$j] + $this->V[$j][$j] * $f;
 
  116                                        for ($k = $j+1; $k <= $i_; ++$k) {
 
  117                                                $g += $this->V[$k][$j] * $this->d[$k];
 
  118                                                $this->
e[$k] += $this->V[$k][$j] * $f;
 
  123                                for ($j = 0; $j < $i; ++$j) {
 
  125                                        $f += $this->
e[$j] * $this->d[$j];
 
  128                                for ($j=0; $j < $i; ++$j) {
 
  129                                        $this->
e[$j] -= $hh * $this->d[$j];
 
  131                                for ($j = 0; $j < $i; ++$j) {
 
  134                                        for ($k = $j; $k <= $i_; ++$k) {
 
  135                                                $this->V[$k][$j] -= ($f * $this->
e[$k] + $g * $this->d[$k]);
 
  137                                        $this->d[$j] = $this->V[$i-1][$j];
 
  138                                        $this->V[$i][$j] = 0.0;
 
  145                for ($i = 0; $i < $this->n-1; ++$i) {
 
  146                        $this->V[$this->n-1][$i] = $this->V[$i][$i];
 
  147                        $this->V[$i][$i] = 1.0;
 
  150                                for ($k = 0; $k <= $i; ++$k) {
 
  151                                        $this->d[$k] = $this->V[$k][$i+1] / 
$h;
 
  153                                for ($j = 0; $j <= $i; ++$j) {
 
  155                                        for ($k = 0; $k <= $i; ++$k) {
 
  156                                                $g += $this->V[$k][$i+1] * $this->V[$k][$j];
 
  158                                        for ($k = 0; $k <= $i; ++$k) {
 
  159                                                $this->V[$k][$j] -= $g * $this->d[$k];
 
  163                        for ($k = 0; $k <= $i; ++$k) {
 
  164                                $this->V[$k][$i+1] = 0.0;
 
  168                $this->d = $this->V[$this->n-1];
 
  169                $this->V[$this->n-1] = array_fill(0, $j, 0.0);
 
  170                $this->V[$this->n-1][$this->n-1] = 1.0;
 
  187                        $this->
e[$i-1] = $this->
e[$i];
 
  189                $this->
e[$this->n-1] = 0.0;
 
  192                $eps  = pow(2.0,-52.0);
 
  196                        $tst1 = max($tst1, abs($this->d[
$l]) + abs($this->
e[$l]));
 
  198                        while ($m < $this->n) {
 
  199                                if (abs($this->
e[$m]) <= $eps * $tst1)
 
  212                                        $p = ($this->d[
$l+1] - $g) / (2.0 * $this->
e[$l]);
 
  216                                        $this->d[
$l] = $this->
e[
$l] / ($p + 
$r);
 
  217                                        $this->d[
$l+1] = $this->
e[
$l] * ($p + 
$r);
 
  218                                        $dl1 = $this->d[
$l+1];
 
  219                                        $h = $g - $this->d[
$l];
 
  227                                        $el1 = $this->
e[$l + 1];
 
  229                                        for ($i = $m-1; $i >= 
$l; --$i) {
 
  233                                                $g  = $c * $this->
e[$i];
 
  235                                                $r  = 
hypo($p, $this->e[$i]);
 
  236                                                $this->e[$i+1] = $s * 
$r;
 
  237                                                $s = $this->e[$i] / 
$r;
 
  239                                                $p = $c * $this->d[$i] - $s * $g;
 
  240                                                $this->d[$i+1] = 
$h + $s * ($c * $g + $s * $this->d[$i]);
 
  243                                                        $h = $this->V[$k][$i+1];
 
  244                                                        $this->V[$k][$i+1] = $s * $this->V[$k][$i] + $c * 
$h;
 
  245                                                        $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * 
$h;
 
  248                                        $p = -$s * $s2 * $c3 * $el1 * $this->
e[
$l] / $dl1;
 
  249                                        $this->
e[
$l] = $s * $p;
 
  250                                        $this->d[
$l] = $c * $p;
 
  252                                } 
while (abs($this->
e[$l]) > $eps * $tst1);
 
  254                        $this->d[
$l] = $this->d[
$l] + $f;
 
  259                for ($i = 0; $i < $this->n - 1; ++$i) {
 
  262                        for ($j = $i+1; $j < 
$this->n; ++$j) {
 
  263                                if ($this->d[$j] < $p) {
 
  269                                $this->d[$k] = $this->d[$i];
 
  272                                        $p = $this->V[$j][$i];
 
  273                                        $this->V[$j][$i] = $this->V[$j][$k];
 
  274                                        $this->V[$j][$k] = $p;
 
  295                for ($m = $low+1; $m <= $high-1; ++$m) {
 
  298                        for ($i = $m; $i <= $high; ++$i) {
 
  299                                $scale = $scale + abs($this->H[$i][$m-1]);
 
  304                                for ($i = $high; $i >= $m; --$i) {
 
  305                                        $this->ort[$i] = $this->H[$i][$m-1] / $scale;
 
  306                                        $h += $this->ort[$i] * $this->ort[$i];
 
  309                                if ($this->ort[$m] > 0) {
 
  312                                $h -= $this->ort[$m] * $g;
 
  313                                $this->ort[$m] -= $g;
 
  316                                for ($j = $m; $j < 
$this->n; ++$j) {
 
  318                                        for ($i = $high; $i >= $m; --$i) {
 
  319                                                $f += $this->ort[$i] * $this->H[$i][$j];
 
  322                                        for ($i = $m; $i <= $high; ++$i) {
 
  323                                                $this->H[$i][$j] -= $f * $this->ort[$i];
 
  326                                for ($i = 0; $i <= $high; ++$i) {
 
  328                                        for ($j = $high; $j >= $m; --$j) {
 
  329                                                $f += $this->ort[$j] * $this->H[$i][$j];
 
  332                                        for ($j = $m; $j <= $high; ++$j) {
 
  333                                                $this->H[$i][$j] -= $f * $this->ort[$j];
 
  336                                $this->ort[$m] = $scale * $this->ort[$m];
 
  337                                $this->H[$m][$m-1] = $scale * $g;
 
  344                                $this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
 
  347                for ($m = $high-1; $m >= $low+1; --$m) {
 
  348                        if ($this->H[$m][$m-1] != 0.0) {
 
  349                                for ($i = $m+1; $i <= $high; ++$i) {
 
  350                                        $this->ort[$i] = $this->H[$i][$m-1];
 
  352                                for ($j = $m; $j <= $high; ++$j) {
 
  354                                        for ($i = $m; $i <= $high; ++$i) {
 
  355                                                $g += $this->ort[$i] * $this->V[$i][$j];
 
  358                                        $g = ($g / $this->ort[$m]) / $this->H[$m][$m-1];
 
  359                                        for ($i = $m; $i <= $high; ++$i) {
 
  360                                                $this->V[$i][$j] += $g * $this->ort[$i];
 
  373        private function cdiv($xr, $xi, $yr, $yi) {
 
  374                if (abs($yr) > abs($yi)) {
 
  377                        $this->cdivr = ($xr + 
$r * $xi) / 
$d;
 
  378                        $this->cdivi = ($xi - 
$r * $xr) / 
$d;
 
  382                        $this->cdivr = (
$r * $xr + $xi) / 
$d;
 
  383                        $this->cdivi = (
$r * $xi - $xr) / 
$d;
 
  404                $eps = pow(2.0, -52.0);
 
  406                $p = $q = 
$r = $s = $z = 0;
 
  410                for ($i = 0; $i < $nn; ++$i) {
 
  411                        if (($i < $low) OR ($i > $high)) {
 
  412                                $this->d[$i] = $this->H[$i][$i];
 
  415                        for ($j = max($i-1, 0); $j < $nn; ++$j) {
 
  416                                $norm = $norm + abs($this->H[$i][$j]);
 
  426                                $s = abs($this->H[
$l-1][
$l-1]) + abs($this->H[
$l][
$l]);
 
  430                                if (abs($this->H[
$l][
$l-1]) < $eps * $s) {
 
  438                                $this->H[
$n][
$n] = $this->H[
$n][
$n] + $exshift;
 
  439                                $this->d[
$n] = $this->H[
$n][
$n];
 
  444                        } 
else if (
$l == 
$n-1) {
 
  446                                $p = ($this->H[
$n-1][
$n-1] - $this->H[
$n][
$n]) / 2.0;
 
  449                                $this->H[
$n][
$n] = $this->H[
$n][
$n] + $exshift;
 
  450                                $this->H[
$n-1][
$n-1] = $this->H[
$n-1][
$n-1] + $exshift;
 
  459                                        $this->d[
$n-1] = 
$x + $z;
 
  460                                        $this->d[
$n] = $this->d[
$n-1];
 
  462                                                $this->d[
$n] = 
$x - 
$w / $z;
 
  464                                        $this->
e[
$n-1] = 0.0;
 
  467                                        $s = abs(
$x) + abs($z);
 
  470                                        $r = sqrt($p * $p + $q * $q);
 
  474                                        for ($j = 
$n-1; $j < $nn; ++$j) {
 
  475                                                $z = $this->H[
$n-1][$j];
 
  476                                                $this->H[
$n-1][$j] = $q * $z + $p * $this->H[
$n][$j];
 
  477                                                $this->H[
$n][$j] = $q * $this->H[
$n][$j] - $p * $z;
 
  480                                        for ($i = 0; $i <= n; ++$i) {
 
  481                                                $z = $this->H[$i][
$n-1];
 
  482                                                $this->H[$i][
$n-1] = $q * $z + $p * $this->H[$i][
$n];
 
  483                                                $this->H[$i][
$n] = $q * $this->H[$i][
$n] - $p * $z;
 
  486                                        for ($i = $low; $i <= $high; ++$i) {
 
  487                                                $z = $this->V[$i][
$n-1];
 
  488                                                $this->V[$i][
$n-1] = $q * $z + $p * $this->V[$i][
$n];
 
  489                                                $this->V[$i][
$n] = $q * $this->V[$i][
$n] - $p * $z;
 
  493                                        $this->d[
$n-1] = 
$x + $p;
 
  494                                        $this->d[
$n] = 
$x + $p;
 
  507                                        $y = $this->H[
$n-1][
$n-1];
 
  513                                        for ($i = $low; $i <= 
$n; ++$i) {
 
  514                                                $this->H[$i][$i] -= 
$x;
 
  516                                        $s = abs($this->H[
$n][
$n-1]) + abs($this->H[
$n-1][
$n-2]);
 
  518                                        $w = -0.4375 * $s * $s;
 
  522                                        $s = (
$y - 
$x) / 2.0;
 
  529                                                $s = 
$x - 
$w / ((
$y - 
$x) / 2.0 + $s);
 
  530                                                for ($i = $low; $i <= 
$n; ++$i) {
 
  531                                                        $this->H[$i][$i] -= $s;
 
  542                                        $z = $this->H[$m][$m];
 
  545                                        $p = (
$r * $s - 
$w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
 
  546                                        $q = $this->H[$m+1][$m+1] - $z - 
$r - $s;
 
  547                                        $r = $this->H[$m+2][$m+1];
 
  548                                        $s = abs($p) + abs($q) + abs(
$r);
 
  555                                        if (abs($this->H[$m][$m-1]) * (abs($q) + abs(
$r)) <
 
  556                                                $eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
 
  561                                for ($i = $m + 2; $i <= 
$n; ++$i) {
 
  562                                        $this->H[$i][$i-2] = 0.0;
 
  564                                                $this->H[$i][$i-3] = 0.0;
 
  568                                for ($k = $m; $k <= 
$n-1; ++$k) {
 
  569                                        $notlast = ($k != 
$n-1);
 
  571                                                $p = $this->H[$k][$k-1];
 
  572                                                $q = $this->H[$k+1][$k-1];
 
  573                                                $r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
 
  574                                                $x = abs($p) + abs($q) + abs(
$r);
 
  584                                        $s = sqrt($p * $p + $q * $q + 
$r * 
$r);
 
  590                                                        $this->H[$k][$k-1] = -$s * 
$x;
 
  591                                                } elseif (
$l != $m) {
 
  592                                                        $this->H[$k][$k-1] = -$this->H[$k][$k-1];
 
  601                                                for ($j = $k; $j < $nn; ++$j) {
 
  602                                                        $p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
 
  604                                                                $p = $p + 
$r * $this->H[$k+2][$j];
 
  605                                                                $this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
 
  607                                                        $this->H[$k][$j] = $this->H[$k][$j] - $p * 
$x;
 
  608                                                        $this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * 
$y;
 
  611                                                for ($i = 0; $i <= min(
$n, $k+3); ++$i) {
 
  612                                                        $p = 
$x * $this->H[$i][$k] + 
$y * $this->H[$i][$k+1];
 
  614                                                                $p = $p + $z * $this->H[$i][$k+2];
 
  615                                                                $this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * 
$r;
 
  617                                                        $this->H[$i][$k] = $this->H[$i][$k] - $p;
 
  618                                                        $this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
 
  621                                                for ($i = $low; $i <= $high; ++$i) {
 
  622                                                        $p = 
$x * $this->V[$i][$k] + 
$y * $this->V[$i][$k+1];
 
  624                                                                $p = $p + $z * $this->V[$i][$k+2];
 
  625                                                                $this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * 
$r;
 
  627                                                        $this->V[$i][$k] = $this->V[$i][$k] - $p;
 
  628                                                        $this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q;
 
  640                for (
$n = $nn-1; 
$n >= 0; --
$n) {
 
  646                                $this->H[
$n][
$n] = 1.0;
 
  647                                for ($i = 
$n-1; $i >= 0; --$i) {
 
  648                                        $w = $this->H[$i][$i] - $p;
 
  650                                        for ($j = 
$l; $j <= 
$n; ++$j) {
 
  651                                                $r = 
$r + $this->H[$i][$j] * $this->H[$j][
$n];
 
  653                                        if ($this->
e[$i] < 0.0) {
 
  658                                                if ($this->
e[$i] == 0.0) {
 
  660                                                                $this->H[$i][
$n] = -
$r / 
$w;
 
  662                                                                $this->H[$i][
$n] = -
$r / ($eps * $norm);
 
  666                                                        $x = $this->H[$i][$i+1];
 
  667                                                        $y = $this->H[$i+1][$i];
 
  668                                                        $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->
e[$i] * $this->
e[$i];
 
  669                                                        $t = (
$x * $s - $z * 
$r) / $q;
 
  670                                                        $this->H[$i][
$n] = 
$t;
 
  671                                                        if (abs(
$x) > abs($z)) {
 
  674                                                                $this->H[$i+1][
$n] = (-$s - 
$y * 
$t) / $z;
 
  678                                                $t = abs($this->H[$i][
$n]);
 
  679                                                if (($eps * 
$t) * 
$t > 1) {
 
  680                                                        for ($j = $i; $j <= 
$n; ++$j) {
 
  681                                                                $this->H[$j][
$n] = $this->H[$j][
$n] / 
$t;
 
  690                                if (abs($this->H[
$n][
$n-1]) > abs($this->H[
$n-1][
$n])) {
 
  691                                        $this->H[
$n-1][
$n-1] = $q / $this->H[
$n][
$n-1];
 
  692                                        $this->H[
$n-1][
$n] = -($this->H[
$n][
$n] - $p) / $this->H[
$n][
$n-1];
 
  694                                        $this->
cdiv(0.0, -$this->H[
$n-1][
$n], $this->H[
$n-1][
$n-1] - $p, $q);
 
  698                                $this->H[
$n][
$n-1] = 0.0;
 
  699                                $this->H[
$n][
$n] = 1.0;
 
  700                                for ($i = 
$n-2; $i >= 0; --$i) {
 
  704                                        for ($j = 
$l; $j <= 
$n; ++$j) {
 
  705                                                $ra = $ra + $this->H[$i][$j] * $this->H[$j][
$n-1];
 
  706                                                $sa = $sa + $this->H[$i][$j] * $this->H[$j][
$n];
 
  708                                        $w = $this->H[$i][$i] - $p;
 
  709                                        if ($this->
e[$i] < 0.0) {
 
  715                                                if ($this->
e[$i] == 0) {
 
  716                                                        $this->
cdiv(-$ra, -$sa, 
$w, $q);
 
  721                                                        $x = $this->H[$i][$i+1];
 
  722                                                        $y = $this->H[$i+1][$i];
 
  723                                                        $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->
e[$i] * $this->
e[$i] - $q * $q;
 
  724                                                        $vi = ($this->d[$i] - $p) * 2.0 * $q;
 
  725                                                        if ($vr == 0.0 & $vi == 0.0) {
 
  726                                                                $vr = $eps * $norm * (abs(
$w) + abs($q) + abs(
$x) + abs(
$y) + abs($z));
 
  728                                                        $this->
cdiv(
$x * 
$r - $z * $ra + $q * $sa, 
$x * $s - $z * $sa - $q * $ra, $vr, $vi);
 
  731                                                        if (abs(
$x) > (abs($z) + abs($q))) {
 
  732                                                                $this->H[$i+1][
$n-1] = (-$ra - 
$w * $this->H[$i][
$n-1] + $q * $this->H[$i][
$n]) / 
$x;
 
  733                                                                $this->H[$i+1][
$n] = (-$sa - 
$w * $this->H[$i][
$n] - $q * $this->H[$i][
$n-1]) / 
$x;
 
  735                                                                $this->
cdiv(-
$r - 
$y * $this->H[$i][
$n-1], -$s - 
$y * $this->H[$i][
$n], $z, $q);
 
  741                                                $t = max(abs($this->H[$i][
$n-1]),abs($this->H[$i][
$n]));
 
  742                                                if (($eps * 
$t) * 
$t > 1) {
 
  743                                                        for ($j = $i; $j <= 
$n; ++$j) {
 
  744                                                                $this->H[$j][
$n-1] = $this->H[$j][
$n-1] / 
$t;
 
  745                                                                $this->H[$j][
$n]   = $this->H[$j][
$n] / 
$t;
 
  754                for ($i = 0; $i < $nn; ++$i) {
 
  755                        if ($i < $low | $i > $high) {
 
  756                                for ($j = $i; $j < $nn; ++$j) {
 
  757                                        $this->V[$i][$j] = $this->H[$i][$j];
 
  763                for ($j = $nn-1; $j >= $low; --$j) {
 
  764                        for ($i = $low; $i <= $high; ++$i) {
 
  766                                for ($k = $low; $k <= min($j,$high); ++$k) {
 
  767                                        $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
 
  769                                $this->V[$i][$j] = $z;
 
  783                $this->A = $Arg->getArray();
 
  784                $this->n = $Arg->getColumnDimension();
 
  801                        $this->ort = array();
 
  817                return new Matrix($this->V, $this->n, $this->n);
 
  851                        $D[$i] = array_fill(0, $this->n, 0.0);
 
  852                        $D[$i][$i] = $this->d[$i];
 
  853                        if ($this->
e[$i] == 0) {
 
  856                        $o = ($this->
e[$i] > 0) ? $i + 1 : $i - 1;
 
  857                        $D[$i][$o] = $this->
e[$i];
 
  859                return new Matrix($D);
 
An exception for terminatinating execution or to throw for unit testing.
getRealEigenvalues()
Return the real parts of the eigenvalues.
cdiv($xr, $xi, $yr, $yi)
Performs complex division.
getV()
Return the eigenvector matrix.
orthes()
Nonsymmetric reduction to Hessenberg form.
hqr2()
Nonsymmetric reduction from Hessenberg to real Schur form.
getImagEigenvalues()
Return the imaginary parts of the eigenvalues.
tred2()
Symmetric Householder reduction to tridiagonal form.
__construct($Arg)
Constructor: Check for symmetry, then construct the eigenvalue decomposition.
tql2()
Symmetric tridiagonal QL algorithm.
getD()
Return the block diagonal eigenvalue matrix.