Subversion Repositories eFlore/Applications.cel

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2388 jpm 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