Subversion Repositories eFlore/Applications.cel

Rev

Details | Last modification | View Log | RSS feed

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