File Coverage

lib/Math/MatrixDecomposition/Eigen.pm
Criterion Covered Total %
statement 530 649 81.6
branch 104 220 47.2
condition 54 120 45.0
subroutine 17 19 89.4
pod 7 8 87.5
total 712 1016 70.0


line stmt bran cond sub pod time code
1             ## Math/MatrixDecomposition/Eigen.pm --- eigenvalues and eigenvectors.
2              
3             # Copyright (C) 2010 Ralph Schleicher. All rights reserved.
4              
5             # This program is free software; you can redistribute it and/or
6             # modify it under the same terms as Perl itself.
7              
8             # Commentary:
9             #
10             # Code derived from EiSPACK procedures and
11             # Jama's 'EigenvalueDecomposition' class.
12              
13             ## Code:
14              
15             package Math::MatrixDecomposition::Eigen;
16              
17 1     1   68965 use strict;
  1         12  
  1         31  
18 1     1   5 use warnings;
  1         1  
  1         30  
19 1     1   4 use Carp;
  1         1  
  1         58  
20 1     1   5 use Exporter qw(import);
  1         2  
  1         34  
21 1     1   5 use POSIX qw(:float_h);
  1         8  
  1         10  
22 1     1   1830 use Math::Complex qw();
  1         12293  
  1         42  
23 1     1   12 use Scalar::Util qw(looks_like_number);
  1         1  
  1         88  
24              
25 1     1   633 use Math::BLAS;
  1         7290  
  1         143  
26 1     1   544 use Math::MatrixDecomposition::Util qw(:all);
  1         2  
  1         181  
27              
28             BEGIN
29             {
30 1     1   5 our $VERSION = '1.06';
31 1         6183 our @EXPORT_OK = qw(eig);
32             }
33              
34             # Calculate eigenvalues and eigenvectors (convenience function).
35             sub eig
36             {
37 0     0 1 0 __PACKAGE__->new (@_);
38             }
39              
40             # Standard constructor.
41             sub new
42             {
43 1     1 1 14 my $class = shift;
44 1         5 my $self =
45             {
46             # Eigenvalues (a vector).
47             value => undef,
48              
49             # Eigenvectors (an array of vectors).
50             vector => undef,
51             };
52              
53 1   33     8 bless ($self, ref ($class) || $class);
54              
55             # Process arguments.
56 1 50       4 $self->decompose (@_)
57             if @_ > 0;
58              
59             # Return object.
60 1         4 $self;
61             }
62              
63             # Calculate eigenvalues and eigenvectors.
64             sub decompose
65             {
66 4     4 1 820 my $self = shift;
67              
68             # Check arguments.
69 4         7 my $a = shift;
70              
71 4 50 33     26 my $m = @_ > 0 && looks_like_number ($_[0]) ? shift : sqrt (@$a);
72 4 50 33     23 my $n = @_ > 0 && looks_like_number ($_[0]) ? shift : $m ? @$a / $m : 0;
    50          
73              
74 4 50 33     79 croak ('Invalid argument')
      33        
      33        
      33        
75             if (@$a != ($m * $n)
76             || mod ($m, 1) != 0 || $m < 1
77             || mod ($n, 1) != 0 || $n < 1);
78              
79 4 50       12 croak ('Matrix has to be square')
80             if $m != $n;
81              
82             # Get properties.
83 4         8 my %prop = @_;
84              
85 4   50     19 $prop{balance} //= 1;
86 4   50     18 $prop{normalize} //= 1;
87 4   50     17 $prop{positive} //= 1;
88              
89             # Index of last row/column.
90 4         8 my $end = $n - 1;
91              
92             # Eigenvalues (a vector).
93 4   100     21 $$self{value} //= [];
94              
95             # Vector $d contains the real part of the eigenvalues and vector $e
96             # contains the imaginary part of the eigenvalues.
97 4         7 my $d = $$self{value};
98 4         7 my $e = [];
99              
100 4 100       12 splice (@$d, $n)
101             if @$d > $n;
102              
103             # Eigenvectors (an array of column vectors).
104 4   100     13 $$self{vector} //= [];
105              
106             # Matrix $Z contains the eigenvectors, please note
107             # that $$Z[$j][$i] denotes matrix element Z(i,j).
108 4         5 my $Z = $$self{vector};
109              
110 4 100       11 splice (@$Z, $n)
111             if @$Z > $n;
112              
113 4         10 for my $v (@$Z)
114             {
115 7 100       17 splice (@$v, $n)
116             if @$v > $n;
117             }
118              
119             # True if matrix is symmetric.
120 4         8 my $sym = 1;
121              
122             SYM:
123              
124 4         12 for my $i (0 .. $end)
125             {
126 6         13 for my $j ($i + 1 .. $end)
127             {
128 6 100       18 if ($$a[$i * $n + $j] != $$a[$j * $n + $i])
129             {
130 3         6 $sym = 0;
131 3         6 last SYM;
132             }
133             }
134             }
135              
136 4 100       11 if ($sym)
137             {
138             # Copy matrix elements.
139 1         2 for my $j (0 .. $end)
140             {
141 3   50     12 $$Z[$j] //= [];
142              
143 3         4 for my $i (0 .. $end)
144             {
145 9         14 $$Z[$j][$i] = $$a[$i * $n + $j];
146             }
147             }
148              
149             # Reduce a real symmetric matrix to a symmetric tridiagonal matrix
150             # using and accumulating orthogonal similarity transformations.
151             #
152             # See EiSPACK procedure 'tred2'.
153 1         2 if (1)
154             {
155 1         2 my ($i, $j, $k, $l,
156             $f, $g, $h, $hh, $scale);
157              
158 1         3 for $i (0 .. $end)
159             {
160 3         4 $$d[$i] = $$Z[$i][$end];
161             }
162              
163 1         15 for $i (reverse (1 .. $end))
164             {
165 2         4 $l = $i - 1;
166 2         2 $h = 0;
167              
168             # Scale row.
169 2         3 $scale = 0;
170              
171 2         4 for $k (0 .. $l)
172             {
173 3         6 $scale += abs ($$d[$k]);
174             }
175              
176 2 50       5 if ($scale == 0)
177             {
178 0         0 $$e[$i] = $$d[$l];
179              
180 0         0 for $j (0 .. $l)
181             {
182 0         0 $$d[$j] = $$Z[$j][$l];
183              
184 0         0 $$Z[$j][$i] = 0;
185 0         0 $$Z[$i][$j] = 0;
186             }
187             }
188             else
189             {
190 2         5 for $k (0 .. $l)
191             {
192 3         5 $$d[$k] /= $scale;
193 3         9 $h += $$d[$k] ** 2;
194             }
195              
196 2         3 $f = $$d[$l];
197 2         8 $g = - sign (sqrt ($h), $f);
198 2         4 $h -= $f * $g;
199              
200 2         4 $$d[$l] = $f - $g;
201 2         3 $$e[$i] = $scale * $g;
202              
203             # Form a*u.
204 2         3 for $j (0 .. $l)
205             {
206 3         5 $$e[$j] = 0;
207             }
208              
209 2         3 for $j (0 .. $l)
210             {
211 3         4 $f = $$d[$j];
212 3         4 $$Z[$i][$j] = $f;
213 3         5 $g = $$e[$j] + $$Z[$j][$j] * $f;
214              
215 3         5 for $k ($j + 1 .. $l)
216             {
217 1         2 $g += $$Z[$j][$k] * $$d[$k];
218 1         2 $$e[$k] += $$Z[$j][$k] * $f;
219             }
220              
221 3         5 $$e[$j] = $g;
222             }
223              
224             # Form p.
225 2         2 $f = 0;
226              
227 2         5 for $j (0 .. $l)
228             {
229 3         4 $$e[$j] /= $h;
230 3         4 $f += $$e[$j] * $$d[$j];
231             }
232              
233             # Form q.
234 2         4 $hh = $f / ($h + $h);
235              
236 2         3 for $j (0 .. $l)
237             {
238 3         13 $$e[$j] -= $hh * $$d[$j];
239             }
240              
241             # Form reduced a.
242 2         3 for $j (0 .. $l)
243             {
244 3         3 $f = $$d[$j];
245 3         5 $g = $$e[$j];
246              
247 3         4 for $k ($j .. $l)
248             {
249 4         8 $$Z[$j][$k] -= ($f * $$e[$k] + $g * $$d[$k]);
250             }
251              
252 3         4 $$d[$j] = $$Z[$j][$l];
253 3         4 $$Z[$j][$i] = 0;
254             }
255             }
256              
257 2         5 $$d[$i] = $h;
258             }
259              
260             # Accumulation of transformation matrices.
261 1         3 for $i (1 .. $end)
262             {
263 2         3 $l = $i - 1;
264              
265 2         3 $$Z[$l][$end] = $$Z[$l][$l];
266 2         3 $$Z[$l][$l] = 1;
267              
268 2         3 $h = $$d[$i];
269 2 50       4 if ($h != 0)
270             {
271 2         4 for $k (0 .. $l)
272             {
273 3         5 $$d[$k] = $$Z[$i][$k] / $h;
274             }
275              
276 2         4 for $j (0 .. $l)
277             {
278 3         3 $g = 0;
279              
280 3         5 for $k (0 .. $l)
281             {
282 5         7 $g += $$Z[$i][$k] * $$Z[$j][$k];
283             }
284              
285 3         4 for $k (0 .. $l)
286             {
287 5         9 $$Z[$j][$k] -= $g * $$d[$k];
288             }
289             }
290             }
291              
292 2         3 for $k (0 .. $l)
293             {
294 3         6 $$Z[$i][$k] = 0;
295             }
296             }
297              
298 1         2 for $j (0 .. $end)
299             {
300 3         4 $$d[$j] = $$Z[$j][$end];
301 3         4 $$Z[$j][$end] = 0;
302             }
303              
304 1         1 $$Z[$end][$end] = 1;
305 1         2 $$e[0] = 0;
306             }
307              
308             # Find the eigenvalues and eigenvectors of a symmetric tridiagonal
309             # matrix by the QL method.
310             #
311             # See EiSPACK procedure 'tql2'.
312 1         2 if (1)
313             {
314 1         3 my ($i, $j, $k, $l, $m,
315             $c, $c2, $c3, $dl1, $el1, $f, $g, $h, $p, $r, $s, $s2, $t, $t2);
316              
317 1         3 for $i (1 .. $end)
318             {
319 2         4 $$e[$i - 1] = $$e[$i];
320             }
321              
322 1         2 $$e[$end] = 0;
323              
324 1         9 $f = 0;
325 1         2 $t = 0;
326              
327 1         3 for $l (0 .. $end)
328             {
329 3         7 $t2 = abs ($$d[$l]) + abs ($$e[$l]);
330 3 100       43 $t = $t2 if $t2 > $t;
331              
332 3         10 for ($m = $l; $m < $n; ++$m)
333             {
334 6 100       12 last if abs ($$e[$m]) <= eps * $t;
335             }
336              
337 3 100       7 if ($m > $l)
338             {
339 2         1 while (1)
340             {
341 5         8 $g = $$d[$l];
342 5         17 $p = ($$d[$l + 1] - $g) / (2 * $$e[$l]);
343 5         10 $r = sign (hypot ($p, 1), $p);
344              
345 5         8 $$d[$l] = $$e[$l] / ($p + $r);
346 5         7 $$d[$l + 1] = $$e[$l] * ($p + $r);
347 5         5 $dl1 = $$d[$l + 1];
348 5         6 $h = $g - $$d[$l];
349              
350 5         8 for $i ($l + 2 .. $end)
351             {
352 4         6 $$d[$i] -= $h;
353             }
354              
355 5         6 $f += $h;
356              
357 5         6 $p = $$d[$m];
358 5         5 $c = 1;
359 5         5 $c2 = $c;
360 5         6 $el1 = $$e[$l + 1];
361 5         5 $s = 0;
362              
363 5         8 for $i (reverse ($l .. $m - 1))
364             {
365 9         10 $c3 = $c2;
366 9         9 $c2 = $c;
367 9         9 $s2 = $s;
368 9         9 $g = $c * $$e[$i];
369 9         9 $h = $c * $p;
370 9         14 $r = hypot ($p, $$e[$i]);
371 9         13 $$e[$i + 1] = $s * $r;
372 9         10 $s = $$e[$i] / $r;
373 9         9 $c = $p / $r;
374 9         10 $p = $c * $$d[$i] - $s * $g;
375 9         11 $$d[$i + 1] = $h + $s * ($c * $g + $s * $$d[$i]);
376              
377 9         11 for $k (0 .. $end)
378             {
379 27         29 $h = $$Z[$i + 1][$k];
380 27         29 $$Z[$i + 1][$k] = $s * $$Z[$i][$k] + $c * $h;
381 27         35 $$Z[$i][$k] = $c * $$Z[$i][$k] - $s * $h;
382             }
383             }
384              
385 5         11 $p = 0 - $s * $s2 * $c3 * $el1 * $$e[$l] / $dl1;
386              
387 5         5 $$e[$l] = $s * $p;
388 5         5 $$d[$l] = $c * $p;
389              
390             # Check convergence.
391 5 100       8 last if abs ($$e[$l]) <= eps * $t;
392             }
393             }
394              
395 3         5 $$d[$l] = $$d[$l] + $f;
396 3         5 $$e[$l] = 0;
397             }
398             }
399             }
400             else
401             {
402             # Hessenberg matrix.
403 3         10 my @H = @$a;
404 3         5 my $N = $n;
405              
406             # Row and column indices of the beginning and end
407             # of the principal sub-matrix.
408 3         6 my $lo = 0;
409 3         5 my $hi = $end;
410              
411             # Permutation vector.
412 3         7 my @perm = ();
413              
414             # Scaling vector.
415 3         5 my @scale = ();
416              
417             # Balance a real matrix and isolate eigenvalues whenever possible.
418 3 50       9 if ($prop{balance})
419             {
420 3         11 my ($p, $s) = gebal ($N, \@H, low => \$lo, high => \$hi);
421              
422 3         6 @perm = @$p;
423 3         8 @scale = @$s;
424             }
425              
426             # Reduce matrix to upper Hessenberg form by orthogonal similarity
427             # transformations.
428             #
429             # See EiSPACK procedure 'orthes'.
430 3         6 my @ort = ();
431              
432 3         22 if (1)
433             {
434 3         6 my ($i, $j, $k, $m,
435             $f, $g, $h, $scale);
436              
437 3         10 for $m ($lo + 1 .. $hi - 1)
438             {
439 2         5 $scale = 0;
440              
441 2         4 for $i ($m .. $hi)
442             {
443 4         9 $scale += abs ($H[$i * $N + ($m - 1)]);
444             }
445              
446 2 50       7 next if $scale == 0;
447              
448 2         2 $h = 0;
449              
450 2         5 for $i (reverse ($m .. $hi))
451             {
452 4         9 $ort[$i] = $H[$i * $N + ($m - 1)] / $scale;
453 4         10 $h += $ort[$i] ** 2;
454             }
455              
456 2         12 $g = - sign (sqrt ($h), $ort[$m]);
457 2         6 $h -= $ort[$m] * $g;
458 2         4 $ort[$m] -= $g;
459              
460 2         4 for $j ($m .. $end)
461             {
462 4         6 $f = 0;
463              
464 4         7 for $i (reverse ($m .. $hi))
465             {
466 8         12 $f += $ort[$i] * $H[$i * $N + $j];
467             }
468              
469 4         4 $f /= $h;
470              
471 4         7 for $i ($m .. $hi)
472             {
473 8         16 $H[$i * $N + $j] -= $f * $ort[$i];
474             }
475             }
476              
477 2         4 for $i (0 .. $hi)
478             {
479 6         9 $f = 0;
480              
481 6         8 for $j (reverse ($m .. $hi))
482             {
483 12         16 $f += $ort[$j] * $H[$i * $N + $j];
484             }
485              
486 6         7 $f /= $h;
487              
488 6         9 for $j ($m .. $hi)
489             {
490 12         16 $H[$i * $N + $j] -= $f * $ort[$j];
491             }
492             }
493              
494 2         4 $ort[$m] *= $scale;
495 2         7 $H[$m * $N + ($m - 1)] = $scale * $g;
496             }
497             }
498              
499             # Accumulate the orthogonal similarity transformations.
500             #
501             # See EiSPACK procedure 'ortran'.
502 3         5 if (1)
503             {
504 3         8 my ($i, $j, $k, $m,
505             $g);
506              
507 3         7 for $j (0 .. $end)
508             {
509 8   100     25 $$Z[$j] //= [];
510              
511 8         14 for $i (0 .. $end)
512             {
513 22 100       37 $$Z[$j][$i] = ($i == $j ? 1 : 0);
514             }
515             }
516              
517 3         15 for $m (reverse ($lo + 1 .. $hi - 1))
518             {
519 2 50       16 next if $H[$m * $N + ($m - 1)] == 0;
520              
521 2         5 for $i ($m + 1 .. $hi)
522             {
523 2         8 $ort[$i] = $H[$i * $N + ($m - 1)];
524             }
525              
526 2         5 for $j ($m .. $hi)
527             {
528 4         6 $g = 0;
529              
530 4         5 for $i ($m .. $hi)
531             {
532 8         12 $g += $ort[$i] * $$Z[$j][$i];
533             }
534              
535 4         10 $g = ($g / $ort[$m]) / $H[$m * $N + ($m - 1)];
536              
537 4         13 for $i ($m .. $hi)
538             {
539 8         24 $$Z[$j][$i] += $g * $ort[$i];
540             }
541             }
542             }
543             }
544              
545             # Find the eigenvalues and eigenvectors of a real upper Hessenberg
546             # matrix by the QR method.
547             #
548             # See EiSPACK procedure 'hqr2'.
549 3         5 if (1)
550             {
551 3         9 my ($i, $j, $k, $l, $m,
552             $h, $g, $f, $p, $q, $r, $s, $t, $w, $x, $y, $z,
553             $norm, $iter, $not_last, $vr, $vi, $ra, $sa);
554              
555 3         6 $t = 0;
556              
557             # Store isolated roots.
558 3         7 for $i (0 .. $end)
559             {
560 8 100 100     25 if ($i < $lo || $i > $hi)
561             {
562 2         5 $$d[$i] = $H[$i * $N + $i];
563 2         3 $$e[$i] = 0;
564             }
565             }
566              
567             # Compute matrix norm.
568 3         6 $norm = 0;
569              
570 3         6 for $i (0 .. $end)
571             {
572 8         12 for $j ($i .. $end)
573             {
574 15         22 $norm += abs ($H[$i * $N + $j]);
575             }
576             }
577              
578             # Search for next eigenvalue.
579 3         4 $iter = 0;
580              
581 3         11 for ($n = $end; $n >= $lo; )
582             {
583             # Look for single small sub-diagonal element.
584 19         37 for ($l = $n; $l > $lo; --$l)
585             {
586 30         40 $s = abs ($H[($l - 1) * $N + ($l - 1)]) + abs ($H[$l * $N + $l]);
587 30 50       42 $s = $norm
588             if $s == 0;
589              
590 30 100       51 last if abs ($H[$l * $N + ($l - 1)]) < eps * $s;
591             }
592              
593 19         31 $x = $H[$n * $N + $n];
594              
595 19 100       36 if ($l == $n)
    100          
596             {
597             # One root found,
598 5         11 $H[$n * $N + $n] = $x + $t;
599              
600 5         24 $$d[$n] = $H[$n * $N + $n];
601 5         9 $$e[$n] = 0;
602              
603 5         5 $n -= 1;
604 5         9 $iter = 0;
605             }
606             elsif ($l == $n - 1)
607             {
608             # Two roots found.
609 1         5 $y = $H[($n - 1) * $N + ($n - 1)];
610 1         4 $w = $H[$n * $N + ($n - 1)] * $H[($n - 1) * $N + $n];
611              
612 1         3 $p = ($y - $x) / 2;
613 1         2 $q = $p * $p + $w;
614 1         3 $z = sqrt (abs ($q));
615              
616 1         3 $H[$n * $N + $n] = $x + $t;
617 1         3 $H[($n - 1) * $N + ($n - 1)] = $y + $t;
618 1         2 $x = $H[$n * $N + $n];
619              
620 1 50       4 if ($q >= 0)
621             {
622             # Real pair.
623 1         34 $z = $p + sign ($z, $p);
624              
625 1         4 $$d[$n - 1] = $x + $z;
626 1         3 $$d[$n] = $$d[$n - 1];
627 1 50       5 $$d[$n] = $x - $w / $z
628             if $z != 0;
629              
630 1         3 $$e[$n - 1] = 0;
631 1         1 $$e[$n] = 0;
632              
633 1         3 $x = $H[$n * $N + ($n - 1)];
634 1         2 $s = abs ($x) + abs ($z);
635 1         2 $p = $x / $s;
636 1         3 $q = $z / $s;
637 1         3 $r = sqrt ($p * $p + $q * $q);
638 1         2 $p = $p / $r;
639 1         1 $q = $q / $r;
640              
641             # Row modification.
642 1         4 for $j ($n - 1 .. $end)
643             {
644 2         4 $z = $H[($n - 1) * $N + $j];
645 2         4 $H[($n - 1) * $N + $j] = $q * $z + $p * $H[$n * $N + $j];
646 2         4 $H[$n * $N + $j] = $q * $H[$n * $N + $j] - $p * $z;
647             }
648              
649             # Column modification.
650 1         3 for $i (0 .. $n)
651             {
652 3         5 $z = $H[$i * $N + ($n - 1)];
653 3         6 $H[$i * $N + ($n - 1)] = $q * $z + $p * $H[$i * $N + $n];
654 3         6 $H[$i * $N + $n] = $q * $H[$i * $N + $n] - $p * $z;
655             }
656              
657             # Accumulate transformations.
658 1         10 for $i ($lo .. $hi)
659             {
660 3         5 $z = $$Z[$n - 1][$i];
661 3         7 $$Z[$n - 1][$i] = $q * $z + $p * $$Z[$n][$i];
662 3         6 $$Z[$n][$i] = $q * $$Z[$n][$i] - $p * $z;
663             }
664             }
665             else
666             {
667             # Complex pair.
668 0         0 $$d[$n - 1] = $x + $p;
669 0         0 $$d[$n] = $x + $p;
670 0         0 $$e[$n - 1] = $z;
671 0         0 $$e[$n] = 0 - $z;
672             }
673              
674 1         2 $n -= 2;
675 1         3 $iter = 0;
676             }
677             else
678             {
679             # Form shift.
680 13         18 $y = $H[($n - 1) * $N + ($n - 1)];
681 13         19 $w = $H[$n * $N + ($n - 1)] * $H[($n - 1) * $N + $n];
682              
683             # Wilkinson's original ad hoc shift.
684 13 50 33     36 if ($iter == 10 || $iter == 20)
685             {
686 0         0 $t += $x;
687              
688 0         0 for $i ($lo .. $n)
689             {
690 0         0 $H[$i * $N + $i] -= $x;
691             }
692              
693 0         0 $s = abs ($H[$n * $N + ($n - 1)]) + abs ($H[($n - 1) * $N + ($n - 2)]);
694 0         0 $x = 0.75 * $s;
695 0         0 $y = $x;
696 0         0 $w = -0.4375 * $s * $s;
697             }
698              
699             # Matlab's new ad hoc shift.
700 13 50       20 if ($iter == 30)
701             {
702 0         0 $s = ($y - $x) / 2;
703 0         0 $s = $s * $s + $w;
704 0 0       0 if ($s > 0)
705             {
706 0         0 $s = sqrt ($s);
707 0 0       0 $s = - $s if $y < $x;
708 0         0 $s = $x - $w / (($y - $x) / 2 + $s);
709              
710 0         0 for $i ($lo .. $n)
711             {
712 0         0 $H[$i * $N + $i] -= $s;
713             }
714              
715 0         0 $t += $s;
716 0         0 $x = 0.964;
717 0         0 $w = $y = $x;
718             }
719             }
720              
721 13         11 ++$iter;
722              
723             # Look for two consecutive small sub-diagonal elements.
724 13         22 for ($m = $n - 2; $m >= $l; --$m)
725             {
726 13         15 $z = $H[$m * $N + $m];
727 13         13 $r = $x - $z;
728 13         12 $s = $y - $z;
729 13         19 $p = ($r * $s - $w) / $H[($m + 1) * $N + $m] + $H[$m * $N + ($m + 1)];
730 13         19 $q = $H[($m + 1) * $N + ($m + 1)] - $z - $r - $s;
731 13         14 $r = $H[($m + 2) * $N + ($m + 1)];
732 13         15 $s = abs ($p) + abs ($q) + abs ($r);
733 13         13 $p = $p / $s;
734 13         12 $q = $q / $s;
735 13         12 $r = $r / $s;
736              
737 13 50       31 last if $m == $l;
738 0 0       0 last if abs ($H[$m * $N + ($m - 1)]) * (abs ($q) + abs ($r)) < eps * (abs ($p) * (abs ($H[($m - 1) * $N + ($m - 1)]) + abs ($z) + abs ($H[($m + 1) * $N + ($m + 1)])));
739             }
740              
741 13         19 for $i ($m + 2 .. $n)
742             {
743 13         17 $H[$i * $N + ($i - 2)] = 0;
744 13 50       26 $H[$i * $N + ($i - 3)] = 0
745             if $i > $m + 2;
746             }
747              
748             # Double QR step.
749 13         18 for $k ($m .. $n - 1)
750             {
751 26         34 $not_last = ($k != $n - 1);
752              
753 26 100       45 if ($k != $m)
754             {
755 13         15 $p = $H[$k * $N + ($k - 1)];
756 13         17 $q = $H[($k + 1) * $N + ($k - 1)];
757 13 50       17 $r = $not_last ? $H[($k + 2) * $N + ($k - 1)] : 0;
758 13         20 $x = abs ($p) + abs ($q) + abs ($r);
759              
760 13 50       18 next if $x == 0;
761              
762 13         14 $p = $p / $x;
763 13         13 $q = $q / $x;
764 13         13 $r = $r / $x;
765             }
766              
767 26         51 $s = sign (sqrt ($p * $p + $q * $q + $r * $r), $p);
768 26 50       37 if ($s != 0)
769             {
770 26 100       43 if ($k != $m)
    50          
771             {
772 13         20 $H[$k * $N + ($k - 1)] = 0 - $s * $x;
773             }
774             elsif ($l != $m)
775             {
776 0         0 $H[$k * $N + ($k - 1)] = - $H[$k * $N + ($k - 1)];
777             }
778              
779 26         371 $p = $p + $s;
780 26         99 $x = $p / $s;
781 26         24 $y = $q / $s;
782 26         24 $z = $r / $s;
783 26         27 $q = $q / $p;
784 26         21 $r = $r / $p;
785              
786 26 100       29 if ($not_last)
787             {
788             # Row modification.
789 13         19 for $j ($k .. $end)
790             {
791 39         58 $p = $H[$k * $N + $j] + $q * $H[($k + 1) * $N + $j] + $r * $H[($k + 2) * $N + $j];
792              
793 39         41 $H[$k * $N + $j] -= $p * $x;
794 39         40 $H[($k + 1) * $N + $j] -= $p * $y;
795 39         47 $H[($k + 2) * $N + $j] -= $p * $z;
796             }
797              
798             # Column modification.
799 13         24 for $i (0 .. min ($n, $k + 3))
800             {
801 39         52 $p = $x * $H[$i * $N + $k] + $y * $H[$i * $N + ($k + 1)] + $z * $H[$i * $N + ($k + 2)];
802              
803 39         41 $H[$i * $N + $k] -= $p;
804 39         43 $H[$i * $N + ($k + 1)] -= $p * $q;
805 39         44 $H[$i * $N + ($k + 2)] -= $p * $r;
806             }
807              
808             # Accumulate transformations.
809 13         20 for $i ($lo .. $hi)
810             {
811 39         56 $p = $x * $$Z[$k][$i] + $y * $$Z[$k + 1][$i] + $z * $$Z[$k + 2][$i];
812              
813 39         40 $$Z[$k][$i] -= $p;
814 39         41 $$Z[$k + 1][$i] -= $p * $q;
815 39         49 $$Z[$k + 2][$i] -= $p * $r;
816             }
817             }
818             else
819             {
820             # Row modification.
821 13         20 for $j ($k .. $end)
822             {
823 26         32 $p = $H[$k * $N + $j] + $q * $H[($k + 1) * $N + $j];
824              
825 26         28 $H[$k * $N + $j] -= $p * $x;
826 26         32 $H[($k + 1) * $N + $j] -= $p * $y;
827             }
828              
829             # Column modification.
830 13         28 for $i (0 .. min ($n, $k + 3))
831             {
832 39         52 $p = $x * $H[$i * $N + $k] + $y * $H[$i * $N + ($k + 1)];
833              
834 39         114 $H[$i * $N + $k] -= $p;
835 39         49 $H[$i * $N + ($k + 1)] -= $p * $q;
836             }
837              
838             # Accumulate transformations.
839 13         18 for $i (0 .. $end)
840             {
841 39         45 $p = $x * $$Z[$k][$i] + $y * $$Z[$k + 1][$i];
842              
843 39         40 $$Z[$k][$i] -= $p;
844 39         57 $$Z[$k + 1][$i] -= $p * $q;
845             }
846             }
847             }
848             }
849             }
850             }
851              
852             # Backsubstitute to find vectors of upper triangular form.
853 3 50       10 return if $norm == 0;
854              
855 3         9 for $n (reverse (0 .. $end))
856             {
857 8         16 $p = $$d[$n];
858 8         10 $q = $$e[$n];
859              
860 8 50       12 if ($q == 0)
    0          
861             {
862             # Real vector.
863 8         8 $m = $n;
864 8         12 $H[$n * $N + $n] = 1;
865              
866 8         15 for $i (reverse (0 .. $n - 1))
867             {
868 7         12 $w = $H[$i * $N + $i] - $p;
869 7         9 $r = 0;
870              
871 7         8 for $j ($m .. $n)
872             {
873 9         17 $r += $H[$i * $N + $j] * $H[$j * $N + $n];
874             }
875              
876 7 50       12 if ($$e[$i] < 0)
877             {
878 0         0 $z = $w;
879 0         0 $s = $r;
880             }
881             else
882             {
883 7         8 $m = $i;
884              
885 7 50       12 if ($$e[$i] == 0)
886             {
887 7 50       15 $H[$i * $N + $n] = ($w != 0 ?
888             0 - $r / $w :
889             0 - $r / (eps * $norm));
890             }
891             else
892             {
893             # Solve real equations.
894 0         0 $x = $H[$i * $N + ($i + 1)];
895 0         0 $y = $H[($i + 1) * $N + $i];
896 0         0 $q = ($$d[$i] - $p) ** 2 + $$e[$i] ** 2;
897 0         0 $t = ($x * $s - $z * $r) / $q;
898              
899 0         0 $H[$i * $N + $n] = $t;
900 0 0       0 $H[($i + 1) * $N + $n] = (abs ($x) > abs ($z) ?
901             (0 - $r - $w * $t) / $x :
902             (0 - $s - $y * $t) / $z);
903             }
904              
905             # Overflow control.
906 7         10 $t = abs ($H[$i * $N + $n]);
907 7 50       14 if ((eps * $t) * $t > 1)
908             {
909 0         0 for $j ($i .. $n)
910             {
911 0         0 $H[$j * $N + $n] /= $t;
912             }
913             }
914             }
915             }
916             }
917             elsif ($q < 0)
918             {
919             # Complex vector.
920 0         0 $m = $n - 1;
921              
922             # Last vector component chosen imaginary so that
923             # eigenvector matrix is triangular.
924 0 0       0 if (abs ($H[$n * $N + ($n - 1)]) > abs ($H[($n - 1) * $N + $n]))
925             {
926 0         0 $H[($n - 1) * $N + ($n - 1)] = $q / $H[$n * $N + ($n - 1)];
927 0         0 $H[($n - 1) * $N + $n] = ($p - $H[$n * $N + $n]) / $H[$n * $N + ($n - 1)];
928             }
929             else
930             {
931 0         0 ($H[($n - 1) * $N + ($n - 1)], $H[($n - 1) * $N + $n])
932             = cdiv (0, - $H[($n - 1) * $N + $n],
933             $H[($n - 1) * $N + ($n - 1)] - $p, $q);
934             }
935              
936 0         0 $H[$n * $N + ($n - 1)] = 0;
937 0         0 $H[$n * $N + $n] = 1;
938              
939 0         0 for $i (reverse (0 .. $n - 2))
940             {
941 0         0 $w = $H[$i * $N + $i] - $p;
942              
943 0         0 $ra = 0;
944 0         0 $sa = 0;
945              
946 0         0 for $j ($m .. $n)
947             {
948 0         0 $ra += $H[$i * $N + $j] * $H[$j * $N + ($n - 1)];
949 0         0 $sa += $H[$i * $N + $j] * $H[$j * $N + $n];
950             }
951              
952 0 0       0 if ($$e[$i] < 0)
953             {
954 0         0 $z = $w;
955 0         0 $r = $ra;
956 0         0 $s = $sa;
957             }
958             else
959             {
960 0         0 $m = $i;
961              
962 0 0       0 if ($$e[$i] == 0)
963             {
964 0         0 ($H[$i * $N + ($n - 1)], $H[$i * $N + $n])
965             = cdiv (- $ra, - $sa, $w, $q);
966             }
967             else
968             {
969             # Solve complex equations.
970 0         0 $x = $H[$i * $N + ($i + 1)];
971 0         0 $y = $H[($i + 1) * $N + $i];
972              
973 0         0 $vr = ($$d[$i] - $p) ** 2 + $$e[$i] ** 2 - $q ** 2;
974 0         0 $vi = ($$d[$i] - $p) * 2 * $q;
975              
976 0 0 0     0 if ($vr == 0 && $vi == 0)
977             {
978 0         0 $vr = eps * $norm * (abs ($w) + abs ($q) + abs ($x) + abs ($y) + abs ($z));
979             }
980              
981 0         0 ($H[$i * $N + ($n - 1)], $H[$i * $N + $n])
982             = cdiv ($x * $r - $z * $ra + $q * $sa,
983             $x * $s - $z * $sa - $q * $ra,
984             $vr,
985             $vi);
986              
987 0 0       0 if (abs ($x) > (abs ($z) + abs ($q)))
988             {
989 0         0 $H[($i + 1) * $N + ($n - 1)] = (0 - $ra - $w * $H[$i * $N + ($n - 1)] + $q * $H[$i * $N + $n]) / $x;
990 0         0 $H[($i + 1) * $N + $n] = (0 - $sa - $w * $H[$i * $N + $n] - $q * $H[$i * $N + ($n - 1)]) / $x;
991             }
992             else
993             {
994 0         0 ($H[($i + 1) * $N + ($n - 1)], $H[($i + 1) * $N + $n])
995             = cdiv (0 - $r - $y * $H[$i * $N + ($n - 1)],
996             0 - $s - $y * $H[$i * $N + $n],
997             $z,
998             $q);
999             }
1000             }
1001              
1002             # Overflow control.
1003 0         0 $t = max (abs ($H[$i * $N + ($n - 1)]), abs ($H[$i * $N + $n]));
1004 0 0       0 if ((eps * $t) * $t > 1)
1005             {
1006 0         0 for $j ($i .. $n)
1007             {
1008 0         0 $H[$j * $N + ($n - 1)] /= $t;
1009 0         0 $H[$j * $N + $n] /= $t;
1010             }
1011             }
1012             }
1013             }
1014             }
1015             }
1016              
1017             # Vectors of isolated roots.
1018 3         9 for $i (0 .. $end)
1019             {
1020 8 100 100     25 if ($i < $lo || $i > $hi)
1021             {
1022 2         3 for $j ($i .. $end)
1023             {
1024 3         7 $$Z[$j][$i] = $H[$i * $N + $j];
1025             }
1026             }
1027             }
1028              
1029             # Multiply by transformation matrix to give
1030             # vectors of original full matrix.
1031 3         8 for $j (reverse ($lo .. $end))
1032             {
1033 7         13 $m = min ($j, $hi);
1034              
1035 7         15 for $i ($lo .. $hi)
1036             {
1037 18         19 $z = 0;
1038              
1039 18         20 for $k ($lo .. $m)
1040             {
1041 36         47 $z += $$Z[$k][$i] * $H[$k * $N + $j];
1042             }
1043              
1044 18         37 $$Z[$j][$i] = $z;
1045             }
1046             }
1047             }
1048              
1049             # Form the eigenvectors of a real general matrix by back
1050             # transforming those of the corresponding balanced matrix
1051             # determined by 'balance'.
1052             #
1053             # See EiSPACK procedure 'balbak'.
1054 3 50       9 if ($prop{balance})
1055             {
1056 3         6 my ($i, $j, $k);
1057              
1058             # Undo scaling.
1059 3         7 for $i ($lo .. $hi)
1060             {
1061 6         8 my $s = $scale[$i];
1062              
1063 6         8 for $j (0 .. $end)
1064             {
1065 18         21 $$Z[$j][$i] *= $s;
1066             }
1067             }
1068              
1069             # Undo permutations.
1070 3         10 for $i (reverse (0 .. $lo - 1))
1071             {
1072 1         2 $k = $perm[$i];
1073 1 50       15 if ($k != $i)
1074             {
1075 0         0 for $j (0 .. $end)
1076             {
1077 0         0 ($$Z[$j][$i], $$Z[$j][$k])
1078             = ($$Z[$j][$k], $$Z[$j][$i]);
1079             }
1080             }
1081             }
1082              
1083 3         11 for $i ($hi + 1 .. $end)
1084             {
1085 1         2 $k = $perm[$i];
1086 1 50       5 if ($k != $i)
1087             {
1088 0         0 for $j (0 .. $end)
1089             {
1090 0         0 ($$Z[$j][$i], $$Z[$j][$k])
1091             = ($$Z[$j][$k], $$Z[$j][$i]);
1092             }
1093             }
1094             }
1095             }
1096             }
1097              
1098             # Create complex eigenvalues.
1099 4         10 for my $i (0 .. $end)
1100             {
1101 11 50       22 $$d[$i] = Math::Complex->make ($$d[$i], $$e[$i])
1102             if $$e[$i] != 0;
1103             }
1104              
1105             # Normalize eigenvectors.
1106             $self->normalize
1107 4 50       29 if $prop{normalize};
1108              
1109             # Make first non-zero vector element a positive number.
1110 4 50       10 if ($prop{positive})
1111             {
1112 4         8 my ($i, $j, $k);
1113              
1114 4         9 for $j (0 .. $end)
1115             {
1116 11         18 for $i (0 .. $end)
1117             {
1118 11 50       20 next if $$Z[$j][$i] == 0;
1119              
1120 11 100       18 if ($$Z[$j][$i] < 0)
1121             {
1122 5         7 for $k ($i .. $end)
1123             {
1124 15         20 $$Z[$j][$k] = - $$Z[$j][$k];
1125             }
1126             }
1127              
1128 11         12 last;
1129             }
1130             }
1131             }
1132              
1133             # Return object.
1134 4         36 $self;
1135             }
1136              
1137             # Balance a real matrix.
1138             #
1139             # See LAPACK procedure 'dgebal'.
1140             sub gebal ($$%)
1141             {
1142             # Arguments.
1143 3     3 0 13 my ($n, $a, %opt) = @_;
1144              
1145             # Default options.
1146 3   50     49 $opt{permute} //= 1;
1147 3   50     14 $opt{scale} //= 1;
1148 3   50     7 $opt{low} //= undef;
1149 3   50     7 $opt{high} //= undef;
1150              
1151             # Index of last row/column.
1152 3         5 my $end = $n - 1;
1153              
1154             # Row and column indices of the beginning and end
1155             # of the principal sub-matrix.
1156 3         11 my $k = 0;
1157 3         6 my $l = $end;
1158              
1159             # Permutation vector.
1160 3         8 my @perm = (0 .. $end);
1161              
1162             # Scaling vector.
1163 3         14 my @scale = map (1, 0 .. $end);
1164              
1165 3 50       11 if ($n > 1)
1166             {
1167             # Isolate eigenvalues.
1168 3 50       7 if ($opt{permute})
1169             {
1170             # Search for rows isolating an eigenvalue
1171             # and push them down.
1172             L:
1173             {
1174 3         4 for my $j (reverse (0 .. $l))
  3         8  
1175             {
1176 7         9 my $swap = 1;
1177              
1178 7         11 for my $i (0 .. $l)
1179             {
1180 20 100 100     50 $swap = 0 if $i != $j && $$a[$j * $n + $i] != 0;
1181             }
1182              
1183 7 100       13 next unless $swap;
1184              
1185 1 50       4 if ($j != $l)
1186             {
1187             # Exchange row and column.
1188 0         0 @perm[$j, $l] = @perm[$l, $j];
1189              
1190 0         0 blas_swap ($l + 1, $a, $a,
1191             x_ind => $j,
1192             x_incr => $n,
1193             y_ind => $l,
1194             y_incr => $n);
1195 0         0 blas_swap ($n - $k, $a, $a,
1196             x_ind => $j * $n + $k,
1197             x_incr => 1,
1198             y_ind => $l * $n + $k,
1199             y_incr => 1);
1200             }
1201              
1202 1         2 $l -= 1;
1203 1         3 next L;
1204             }
1205             }
1206              
1207             # Search for columns isolating an eigenvalue
1208             # and push them left.
1209             K:
1210             {
1211 3         6 for my $j ($k .. $l)
  3         9  
1212             {
1213 7         9 my $swap = 1;
1214              
1215 7         11 for my $i ($k .. $l)
1216             {
1217 19 100 66     51 $swap = 0 if $i != $j && $$a[$i * $n + $j] != 0;
1218             }
1219              
1220 7 100       16 next unless $swap;
1221              
1222 1 50       2 if ($j != $k)
1223             {
1224             # Exchange row and column.
1225 0         0 @perm[$j, $k] = @perm[$k, $j];
1226              
1227 0         0 blas_swap ($l + 1, $a, $a,
1228             x_ind => $j,
1229             x_incr => $n,
1230             y_ind => $k,
1231             y_incr => $n);
1232 0         0 blas_swap ($n - $k, $a, $a,
1233             x_ind => $j * $n + $k,
1234             x_incr => 1,
1235             y_ind => $k * $n + $k,
1236             y_incr => 1);
1237             }
1238              
1239 1         3 $k += 1;
1240 1         2 next K;
1241             }
1242             }
1243             }
1244              
1245             # Balance the sub-matrix in rows k to l.
1246 3 50       9 if ($opt{scale})
1247             {
1248 3         8 my ($b, $c, $ca, $f, $g, $r, $ra, $s, $no_conv,
1249             $sfmin1, $sfmax1, $sfmin2, $sfmax2);
1250              
1251             # Scale factors are powers of two.
1252 3         4 $b = 2;
1253              
1254 3         5 if (1)
1255             {
1256 3         10 my $base = 2;
1257 3         6 my $eps = DBL_EPSILON;
1258 3         5 my $small = 1 / DBL_MAX;
1259             # Safe minimum, such that 1/sfmin does not overflow.
1260 3 50       7 my $sfmin = ($small >= DBL_MIN ? $small * (1 + $eps) : DBL_MIN);
1261              
1262 3         5 $sfmin1 = $sfmin / ($eps * $base);
1263 3         5 $sfmax1 = 1 / $sfmin1;
1264 3         4 $sfmin2 = $sfmin1 * $b;
1265 3         5 $sfmax2 = 1 / $sfmin2;
1266             }
1267              
1268             # Iterative loop for norm reduction.
1269 3         3 while (1)
1270             {
1271 4         5 $no_conv = 0;
1272              
1273 4         10 for my $i ($k .. $l)
1274             {
1275 9         27 $c = blas_norm ($l - $k + 1, $a,
1276             norm => BLAS_TWO_NORM,
1277             x_ind => $k * $n + $i,
1278             x_incr => $n);
1279 9         383 $r = blas_norm ($l - $k + 1, $a,
1280             norm => BLAS_TWO_NORM,
1281             x_ind => $i * $n + $k,
1282             x_incr => 1);
1283 9         327 (undef, $ca) = blas_amax_val ($l + 1, $a,
1284             x_ind => $i,
1285             x_incr => $n);
1286 9         212 (undef, $ra) = blas_amax_val ($n - $k, $a,
1287             x_ind => $i * $n + $k,
1288             x_incr => 1);
1289              
1290             # Guard against zero c or r due to underflow.
1291 9 50 33     261 next if $c == 0 || $r == 0;
1292              
1293 9         13 $s = $c + $r;
1294 9         9 $f = 1;
1295              
1296 9         11 $g = $r / $b;
1297 9   66     27 while ($c < $g
      66        
1298             && max ($f, max ($c, $ca)) < $sfmax2
1299             && min ($g, min ($r, $ra)) > $sfmin2)
1300             {
1301 1 50       6 croak ('Infinite loop')
1302             if isnan ($f + $c + $ca + $g + $r + $ra);
1303              
1304 1         2 $f *= $b;
1305 1         4 $c *= $b;
1306 1         2 $ca *= $b;
1307 1         2 $g /= $b;
1308 1         2 $r /= $b;
1309 1         3 $ra /= $b;
1310             }
1311              
1312 9         10 $g = $c / $b;
1313 9   66     20 while ($g >= $r
      66        
1314             && max ($r, $ra) < $sfmax2
1315             && min ($g, min ($f, min ($c, $ca))) > $sfmin2)
1316             {
1317 2         4 $f /= $b;
1318 2         3 $c /= $b;
1319 2         2 $ca /= $b;
1320 2         3 $g /= $b;
1321 2         2 $r *= $b;
1322 2         4 $ra *= $b;
1323             }
1324              
1325             # Now balance.
1326 9 50 66     65 if (! (($c + $r) >= 0.95 * $s
      33        
      66        
      66        
      33        
      33        
1327             || ($f < 1 && $scale[$i] < 1 && $f * $scale[$i] <= $sfmin1)
1328             || ($f > 1 && $scale[$i] > 1 && $scale[$i] >= $sfmax1 / $f)))
1329             {
1330 3         12 blas_rscale ($n - $k, $a,
1331             alpha => $f,
1332             x_ind => $i * $n + $k,
1333             x_incr => 1);
1334 3         96 blas_scale ($l + 1, $a,
1335             alpha => $f,
1336             x_ind => $i,
1337             x_incr => $n);
1338              
1339 3         45 $scale[$i] *= $f;
1340 3         6 $no_conv = 1;
1341             }
1342             }
1343              
1344 4 100       11 last unless $no_conv;
1345             }
1346             }
1347             }
1348              
1349 3         9 ${$opt{low}} = $k
1350 3 50       12 if ref ($opt{low});
1351              
1352 3         4 ${$opt{high}} = $l
1353 3 50       11 if ref ($opt{high});
1354              
1355 3         11 return (\@perm, \@scale);
1356             }
1357              
1358             # Normalize eigenvectors.
1359             sub normalize
1360             {
1361 4     4 1 9 my $self = shift;
1362              
1363             # Work variables.
1364 4         6 my $len;
1365              
1366 4         6 for my $v (@{ $$self{vector} })
  4         14  
1367             {
1368 11         257 $len = blas_norm (@$v, $v, norm => BLAS_TWO_NORM);
1369 11 50       507 blas_rscale (@$v, $v, alpha => $len) if $len != 0;
1370             }
1371              
1372             # Return object.
1373 4         106 $self;
1374             }
1375              
1376             # Sort eigenvalues and corresponding eigenvectors.
1377             sub sort
1378             {
1379 4     4 1 8 my $self = shift;
1380 4   50     14 my $order = shift // 'abs_desc';
1381              
1382             # Permutation vector.
1383 4         6 my @p = ();
1384              
1385 4 50       22 if ($order =~ m/\Avec_/)
    50          
1386             {
1387 0         0 my $Z = $$self{vector};
1388              
1389             @p = ($order eq 'vec_desc' ?
1390 0         0 sort { _cmp_vec ($$Z[$b], $$Z[$a]) } 0 .. $#$Z :
1391             ($order eq 'vec_asc' ?
1392 0 0       0 sort { _cmp_vec ($$Z[$a], $$Z[$b]) } 0 .. $#$Z :
  0 0       0  
1393             croak ("Invalid argument")));
1394             }
1395 4         27 elsif (grep (ref ($_), @{ $$self{value} }))
1396             {
1397             # Consider complex eigenvalues.
1398 0         0 my (@d, @e, @m) = ();
1399              
1400 0         0 for (@{ $$self{value} })
  0         0  
1401             {
1402 0 0       0 if (ref ($_))
1403             {
1404 0         0 push (@d, $_->Re);
1405 0         0 push (@e, $_->Im);
1406 0         0 push (@m, abs ($_));
1407             }
1408             else
1409             {
1410 0         0 push (@d, $_);
1411 0         0 push (@e, 0);
1412 0         0 push (@m, abs ($_));
1413             }
1414             }
1415              
1416             @p = ($order eq 'abs_desc' ?
1417 0 0 0     0 sort { abs ($d[$b]) <=> abs ($d[$a]) || $m[$b] <=> $m[$a] || $d[$b] <=> $d[$a] || $e[$b] <=> $e[$a] } 0 .. $#d :
      0        
1418             ($order eq 'abs_asc' ?
1419 0 0 0     0 sort { abs ($d[$a]) <=> abs ($d[$b]) || $m[$a] <=> $m[$b] || $d[$a] <=> $d[$b] || $e[$a] <=> $e[$b] } 0 .. $#d :
      0        
1420             ($order eq 'norm_desc' ?
1421 0 0 0     0 sort { $m[$b] <=> $m[$a] || $d[$b] <=> $d[$a] || $e[$b] <=> $e[$a] } 0 .. $#d :
1422             ($order eq 'norm_asc' ?
1423 0 0 0     0 sort { $m[$a] <=> $m[$b] || $d[$a] <=> $d[$b] || $e[$a] <=> $e[$b] } 0 .. $#d :
1424             ($order eq 'desc' ?
1425 0 0       0 sort { $d[$b] <=> $d[$a] || $e[$b] <=> $e[$a] } 0 .. $#d :
1426             ($order eq 'asc' ?
1427 0 0       0 sort { $d[$a] <=> $d[$b] || $e[$a] <=> $e[$b] } 0 .. $#d :
  0 0       0  
    0          
    0          
    0          
    0          
    0          
1428             croak ("Invalid argument")))))));
1429             }
1430             else
1431             {
1432             # Only real eigenvalues.
1433 4         10 my $d = $$self{value};
1434              
1435             @p = ($order eq 'abs_desc' || $order eq 'norm_desc' ?
1436 0 0       0 sort { abs ($$d[$b]) <=> abs ($$d[$a]) || $$d[$b] <=> $$d[$a] } 0 .. $#$d :
1437             ($order eq 'abs_asc' || $order eq 'norm_asc' ?
1438 0 0       0 sort { abs ($$d[$a]) <=> abs ($$d[$b]) || $$d[$a] <=> $$d[$b] } 0 .. $#$d :
1439             ($order eq 'desc' ?
1440 0         0 sort { $$d[$b] <=> $$d[$a] } 0 .. $#$d :
1441             ($order eq 'asc' ?
1442 4 50 33     64 sort { $$d[$a] <=> $$d[$b] } 0 .. $#$d :
  9 50 33     37  
    50          
    50          
1443             croak ("Invalid argument")))));
1444             }
1445              
1446             # Reorder eigenvalues and corresponding eigenvectors.
1447 4 50       27 if (@p > 0)
1448             {
1449 4         6 my $ref;
1450              
1451 4         8 $ref = $$self{value};
1452 4         19 @$ref = @$ref[@p];
1453              
1454 4         8 $ref = $$self{vector};
1455 4         18 @$ref = @$ref[@p];
1456             }
1457              
1458             # Return object.
1459 4         10 $self;
1460             }
1461              
1462             # Compare two eigenvectors.
1463             sub _cmp_vec
1464             {
1465 0     0   0 my ($u, $v) = @_;
1466              
1467 0 0       0 croak ("Invalid argument")
1468             if $#$u != $#$v;
1469              
1470 0         0 my $d = 0;
1471              
1472 0         0 for my $i (0 .. $#$u)
1473             {
1474 0 0       0 last if $d = ($$u[$i] <=> $$v[$i]);
1475             }
1476              
1477 0         0 $d;
1478             }
1479              
1480             # Return one or more eigenvalues.
1481             sub value
1482             {
1483 4     4 1 9 my $self = shift;
1484              
1485 4 50       12 @_ > 0 ? @{ $$self{value} }[@_] : @{ $$self{value} };
  0         0  
  4         15  
1486             }
1487              
1488             *values = \&value;
1489              
1490             # Return one or more eigenvectors.
1491             sub vector
1492             {
1493 11     11 1 3797 my $self = shift;
1494              
1495 11 50       30 @_ > 0 ? @{ $$self{vector} }[@_] : @{ $$self{vector} };
  11         32  
  0            
1496             }
1497              
1498             *vectors = \&vector;
1499              
1500             1;
1501              
1502             __END__