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